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Abstract 


A unified construction of high order shape functions is given for all four classical energy 
spaces {H^, H (curl), H (div) and and for elements of “all” shapes (segment, quadrilateral, 
triangle, hexahedron, tetrahedron, triangular prism and pyramid). The discrete spaces spanned 
by the shape functions satisfy the commuting exact sequence property for each element. The 
shape functions are conforming, hierarchical and compatible with other neighboring elements 
across shared boundaries so they may be used in hybrid meshes. Expressions for the shape 
functions are given in coordinate free format in terms of the relevant affine coordinates of each 
element shape. The polynomial order is allowed to differ for each separate topological entity 
(vertex, edge, face or interior) in the mesh, so the shape functions can be used to implement local 
p adaptive finite element methods. Each topological entity may have its own orientation, and 
the shape functions can have that orientation embedded by a simple permutation of arguments. 
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1 Preliminaries 


1.1 Introduction 

In the context of finite elements, construction of higher order shape functions for elements forming the 
exact sequence has been a long standing activity in both engineering and numerical analysis communities. 
A comprehensive review of the subject can be found for example in Demkowicz (2006), Demkowicz et al. 
(2007) and references therein. 

This document presents a self-contained systematic theory for the construction of a particular set of 
hierarchical, orientation embedded, 77(curl), 77(div), and conforming shape functions for elements 

of “all shapes”, forming the ID, 2D, and 3D commuting exact sequences discussed within. By elements of 
“all shapes”, we specifically mean the segment (unit interval) in ID, the quadrilateral and triangle in 2D, 
and the hexahedron, tetrahedron, prism (wedge) and pyramid in 3D. 



Segment Quadrilateral Triangle Hexahedron Tetrahedron Prism Pyramid 

Figure 1.1: Elements of “all shapes”. 


There are many ways to construct sets of shape functions satisfying the aforementioned properties. 
However, we believe that in this work we have constructed a set which strikes an uncommon balance be¬ 
tween simplicity and applicability. For all elements, and each associated energy space, we rely upon a 
simple methodology and a very small collection of ancillary functions to generate all of our shape functions. 
Furthermore, we have supplemented this text with a package written in Fortran 90 defining each function 
presented in this work.^ For these reasons, when reproducing our work in their own software, the readers 
should find the burden of implementation minimal. We hope that our exposition will be clear and useful, 
particularly to those less familiar with the subject. 

For those at the forefront of shape function construction, we hope that our work will be intriguing if only 
for the elegance of our construction. Particularly, we evidence §9 on pyramid shape functions. The higher 
order discrete commuting exact sequence for this element appeared only recently in the work of Nigam and 
Phillips (2012). Our construction for the pyramid presents shape functions spanning each of their discrete 
energy spaces while maintaining compatibility with the other 3D elements. We also remark that, for any 
given mesh, our shape functions are fully compatible across adjacent interelement boundaries due to con¬ 
sidering so-called orientation embeddings. Hence, no alterations of the shape functions are necessary at the 
finite element assembly procedure. Moreover, these orientation embeddings are handled almost effortlessly 
by simply permuting the entries of a few relevant functions. 

With regard to the choice of geometry (shape and size) of master elements, we followed Demkowicz 
*See the ESEAS library available at https://github.com/libESEAS/ESEAS. 
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(2006). However, one of the key points is that our eonstruetion naturally applies to any other ehoiee of 
master element geometries. Other speeifie ehoiees we made when enumerating element vertiees, edges, and 
faees, ean also be modified with little effort to the preferenees of the reader. 

For eompleteness, we have ehosen to thoroughly verify the mathematieal properties and to give a sound 
geometrieal interpretation of our eonstruetions rather than only give the neeessary shape funetions to the 
reader. We eoneede that due to the depth of our presentation, and the expanse of our eoverage, our offense 
lies only in the length of this doeument. However, in a sense, an abridged version of this work is already 
present in a set of tables summarizing all the shape funetions. These ean be eonveniently eonsulted in 
Appendix E. 

1.2 Energy Spaces and Exact Sequences 

Let n C with = 1, 2,3 be a domain. One arrives naturally at the energy spaces iF^(O), iF(eurl, O), 
if(div, O) and L^(0) in eontext of various variational formulations, see e.g. Chapter 1 in Demkowiez et al. 
(2007), Demkowiez and Gopalakrishnan (2014) and Demkowiez (2015). Along with operations of gradient, 
curl and divergence (understood in the sense of distributions), these spaces form the so-called complexes, 
i.e. the composition of any two operators in the sequence reduces to the trivial operator. The ID complex, 
where C M, provides the simplest example: 

M A A l2(S2) a {0} . 

Here, the symbol M denotes constant functions, and {0} is the trivial vector space consisting of the zero 
function only. By using the name of complex, we communicate two simple facts: a) the derivative of a 
constant function is zero, and b) the composition of derivative (in fact, any linear operator) with the trivial 
(zero) operator is trivial as well. Equivalently, we can express the same facts by using null spaces and ranges 
of the involved operators: 

R(id) C N(V) and R(V) C N(0). 

If instead of inclusions above, we have equalities, then we say that the complex (sequence) is exact. This 
is indeed the case for the simply connected domain = (0,1). By using the name exact sequence, we 
communicate more information: a) the derivative of a function is zero if and only if the function is a 
constant, and b) the function V : —)■ is a surjection (onto). Erom now on, we remove mention 

of the first and final terms of fhe exacf sequence. The fwo spaces M and {0}, and fhe operators id and 0, are 
always assumed fo buffress each of fhe sequences we later presenf. Moreover, whenever possible, we absorb 
fhe n assignmenf wifhin fhe nofafion of each energy space. The domain will always be assumed fo be a 
simply connecfed domain in fhe relevanf 

ID Exact Sequence. We now presenf fhe firsf exacf sequence of simply connecfed domains in M: 

( 1 . 1 ) 
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2D Exact Sequence. The exact sequence for simply connected domains in is of the form 

^ ^ , ( 1 . 2 ) 

where V x and x are understood in two dimensions: 


V X E = \7 X 



dE2 _ dEi 
d^i d^2 


ExE 



E 1 E 2 — E 2 F 1 . (1-3) 


By “rotating” i?(curl), the space i7(div) arises naturally: 

ff(div) = [Ve = {M)e F(curl)} . (1.4) 


Defined in this way, the “rotated” exact sequence is immediately satisfied: 

Ri cu^ H{d\v) ^ , (1.5) 

where, for all cj) £ and all E G Tf(curl), the operations satisfy the following relations: 

0 1\ / 0 A 

]V(p, V-Ve = V-[ ]e = VxE. (1.6) 

-1 oj \-l oj 

3D Exact Sequence. Finally, for a simply connected domain in M^, we have the 3D exact sequence 

^ i7(curl) ^ H{dw) ^ . (1.7) 


curl {(!)) = 



For all elements, these exact sequences will be reproduced on the discrete level by replacing the energy 
spaces with appropriate polynomial subspaces.^ We shall use the standard notation: 


ID : 

WP 

V 

- ^ 

yp 




f WP 

V 

—> 

QP 

Vx 
—)■ 

yp 

2D : 






1 WP 

curl 
—> 

yp 

V- 

—^ 

yp 

3D : 

WP 

V 

—^ 

QP 

Vx 
-)■ 

yp 


(1.8) 


The symbol p loosely denotes the polynomial order and should not be interpreted literally.^ 

In this work, we shall consider only spaces of the first type which, with the exception of the pyramid, 
were introduced by Nedelec (1980) in the^zr^t of his two famous papers. The pyramid spaces were taken as 
the first set of spaces proposed by Nigam and Phillips (2012). All these spaces satisfy a number of funda¬ 
mental properties. First, the spaces of the different elements are said be tracewise compatible at the level of 

^Or rational polynomial subspaces in the case of the pyramid (see §9). 

^By this we mean that p should, in fact, be interpreted as a multi-index for the Cartesian product elements. 
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spaces. This allows them to be used in hybrid meshes, which may contain elements of all shapes. Second, 
for each element, contains polynomials of total order p, while Q^, and contain polynomials of 
total order p — 1, meaning that the overall drop in polynomial degree from the first discrete energy space in 
the exact sequence to the last discrete energy space is one. Thirdly, for a given element and energy space, 
the discrete spaces form a nested sequence of spaces as the order increases (for instance, C and 

so on). This is a necessary condition for the construction of hierarchical sets of shape functions. Lastly, 
the spaces form commuting exact sequences for each element. This, coupled with the previous properties, 
ensures (global) interpolation estimates for all of our energy spaces subject to affine transformations of the 
master element geometries (Monk and Demkowicz, 2000). 

1.3 Shape Functions 

We shall always identify the discrete spaces first (like W^, Qp, and in (1.8)), and only afterward 
introduce the corresponding shape functions that provide bases for those spaces. This is a good place to 
remind the reader that there are, in fact, two competing schools of thought when it comes to the theory of 
shape functions. 

The classical definition of Ciarlet (1994) starts with degrees of freedom that are functionals defined on 
some large subset X of an energy space U (like n C H^). The shape functions, which are elements 
spanning some discrete (finite dimensional) space X Y X (e.g. C C°^ H H^), are then defined as 
the dual basis to the linearly independent (when restricted to X) degrees of freedom. An interpolation 
operator from Tf to A is then naturally defined. In this construction, we must (usually) precompute the 
shape functions, e.g. in terms of combinations of monomials whose corresponding coefficients are stored. 

The competing approach of Szabo (Szabo and Babuska, 1991) starts with a direct construction of shape 
functions by following a topological classification (of vertices, edges, faces and element interiors) induced 
by conformity requirements. The shape functions are defined in terms of families of polynomials (e.g. 
Legendre) and their integrals, and are computed using simple recursive formulas. This is the approach taken 
in this work. The so-called projection-based interpolation (Demkowicz, 2006; Demkowicz et ah, 2007) 
defined through local projections over element edges, faces and interior, is introduced independently of the 
construction of shape functions. Therefore, with no need to precompute coefficients defining the shape 
functions, following Szabo’s approach is perhaps more convenient and straightforward. 

1.4 Hierarchy in p 

Given an energy space (like and a conforming discrete space of order p (like C H^), denote the 
shape functions forming a basis for the discrete space by (e.g. span(.BP) = C A construction 
is said to be hierarchical in p if B^ C B^^^ for all p, so that the set of shape functions spanning a space of a 
certain order is found in all subsequent enriched spaces of higher order. This implies that as p increases, all 
one has to do is to add a few functions to a smaller previously constructed set of shape functions. 
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In our construction, hierarchy in p will be enforced. In a given mesh, it will allow comparison of shape 
functions between adjacent elements that have different order, so that at least some of the shape functions of 
the neighboring elements will match. This is crucial with regard to the notion of local p adaptivity, which 
some methods employ. In our work, this flexibility in the variability of the order will be partly reflected by 
the natural anisotropies present in Cartesian product elements, such as quadrilaterals, hexahedra and prisms, 
where each independent direction can have a different order. More information can be consulted in the 
literature (see Demkowicz et al. (2007) and references therein). 

1.5 Traces and Compatibility 

For a function to be contained in a given energy space it must satisfy some global conformity conditions 
which depend upon the space. For instance, functions in are almost globally continuous, but functions 
in can be much more discontinuous. 

Due to these conformity requirements, each energy space has a different definition of trace at the bound¬ 
aries. The different traces only make sense on certain parts of the boundary. For instance, consider a poly¬ 
hedral element. 

• The trace is the value of the function itself at the boundary. In 3D, it may take values at element 
vertices, edges and faces which lie along the boundary. 

• The H (curl) (tangential) trace is the tangential component of the vector valued function across the 
boundary. It may take values at edges and faces in the boundary, but not at vertices, since these do not 
have a concept of tangent. In fact, the H (curl) (tangential) trace is scalar valued across edges (which 
have ID tangent spaces) and has two components across faces (which have 2D tangent spaces). 

• The H (div) (normal) trace is the normal component of the vector valued function across the boundary. 
In 3D, it can take values at the faces of the boundary, but not at vertices or edges, since they do not 
have a unique notion of normal. In 2D, edges along the boundary do have a notion of normal, so the 
(normal) edge trace does exist. 

• There is no notion of trace for L?. 

At the discrete level, these considerations lead naturally to a classification of shape functions according 
to topological entities, which in turn depend on the number of spatial dimensions: 

ID: vertex and edge shape functions, 

2D: vertex, edge, and face shape functions, 

3D: vertex, edge, face and interior shape functions. 

For any given mesh and energy space, shape functions of adjacent elements must be continuous at the 
trace level across the shared interelement boundaries. This is referred to as compatibility, and results in the 
global conformity of the (disjoint union of) shape functions. For instance, in 2D, the order p edge shape 
function of a quadrilateral would need to be compatible with the order p edge shape function of an adjacent 
triangle (see Figure 1.2). 
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Quadrilateral 
edge shape function 

Compatibility: ^ 


Triangle 

edge shape function 



Global HI function 

Global conformity: 


Figure 1.2: Example of compatible edge functions resulting in a globally conforming function. 

1.6 Embedded Sequences and Dimensional Hierarchy 


To enforce compatibility, it is useful to actually begin with a known trace over the (shared) boundary and 
then extend (or lift) it to the rest of the element. This is the approach inherently present in our constructions. 
In fact, this idea is reinforced when looking at the exact sequences. Note the crucial fact that the lower 
dimensional sequences are “embedded” in the higher dimensional sequences if one considers the appropriate 
restrictions. This is better represented by the following diagram. 


3D: 

2D: 

ID: 


^ if (curl) ii(div) L2 


tr 

1 


tr 


^ ii(curl) 

tr 


Vx 




tr 

2 


H 


tr 

1 V 


V 


(1.9) 


tr 

where the “mapping arrows”, , indicate that the range (of the trace) is actually a larger space. These 

arrows are meant to be motivational only. At the discrete level, we always reproduce the above diagram 
precisely with the dotted lines being replaced by well defined maps. 

Indeed, we will enforce a dimensional hierarchy through traces which is consistent with the previous 
discussion. The template for this new form of hierarchy is precisely (1.9), but it is satisfied at the level of the 
shape functions themselves (not only the spaces). One will begin the construction by first defining the ID 
shape functions (over the segment), then defining all the 2D functions (over the quadrilateral and triangle), 
and finish with the 3D shape functions. Throughout the construction, the higher dimensional shape functions 
will have as (nonzero) trace a lower dimensional function, so that they are actually extensions of these lower 
dimensional functions. Indeed, the shape functions are nested through the trace operation at each topological 
entity lying in the boundary. 

For example, given an edge of a 2D element, the nonzero edge traces of the 2D shape functions 
should reproduce the ID shape functions. Hence, some 2D shape functions are said to be extensions 
of the ID shape functions. Similarly, the nonzero edge trace of 2D iF(curl) shape functions should be 
the ID shape functions (see (1.9)). These relations hold per topological entity as well. For example, the 
nonzero edge trace of 2D vertex functions should coincide with the ID vertex functions, and the 
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nonzero trace of the 2D edge functions should coincide with the ID edge functions (see Figures 1.5 
and 1.6 later on). When enforced, these nice relationships not only aid in the compatibility, but, from the 
computational standpoint, have the benefit of allowing us to recycle a large amount of code when moving 
from one element construction to another. 

1.7 Basic Properties of Shape Functions 

In this section we will describe the basic properties that our shape functions should satisfy. These properties 
are specific fo fhe topological enfify (vertices, edges, faces, inferior) and fhe energy space fo which fhe shape 
funcfions are associafed. For example, verfex funclions satisfy differenl properties fhan edge functions in 
and 2D edge functions satisfy differenl properties in and Ff(curl). 

For a given elemenl and energy space, each lopological entity owns a sef of shape functions, which is 
said fo be associated fo fhe enfify. These funclions are nonzero al fhe associafed topological enfify. Indeed, 
each verfex is associafed fo one verfex shape funclion. Meanwhile, each edge, face, and inferior of fhe 
elemenl is associafed fo a set of edge, face, and inferior shape functions of size of fhe order of p, p^ and p^ 
respectively. For example, in each edge is associafed fo a sef of p — 1 edge shape functions. 

Now, for a given dimension, space, and lopological enfify, we will cover Ihree main aspecls of fhe shape 
funclions. Firsl, are fhe vanishing properlies lhal Ihey should satisfy. These properlies eslablish a form of 
Irivial compalibilily along some parls of fhe boundary. Funclions whose Irace vanishes everywhere along fhe 
boundary are called bubbles, and Ihey are Irivially compatible wilh each olher. Second, come fhe nonzero 
Irace properlies. These are, in general, nonlrivial, and ensure eilher full compalibilily or compalibilily 
modulo “orienfalions” (see §1.8 for a discussion on “orienfalions”). Forlunalely, dimensional hierarchy will 
delermine fhe form of Ihese nonzero Iraces. Third, are some properlies lhal fhe shape funclions should 
salisfy along fhe elemenl ilself in order fo have hierarchy in p. 

1.7.1 ID 

In one dimension, classification of shape funclions is: 

: vertex and edge functions, 

; edge functions. 

There is only one simply connected ID elemenl, which is fhe segmenl (or edge), and ils boundary is jusl ils 
Iwo verlices. 

H^. Vertex functions: Firsl, Ihere are fhe vanishing properties. The verfex funclions should vanish al Ihe 
olher (unassocialed) vertex. Second, Ihere are Ihe nonzero Irace properties. The vertex functions should lake 
Ihe value 1 al Ihe associated vertex. This will ensure full compatibility in ID. Third, Ihere is Ihe form of 
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the function itself. There are multiple ways in which the vertex function can decay towards the other vertex 
(see Figure 1.3). If the hierarchy in p is not an issue, the decay could be nonlinear and dependent on p, 
and this may have some computational advantages. For example, when having a mesh with uniform order 
p across all elements, this nonlinear decay might lead to a better conditioning of finite element matrices. 
However, as mentioned before, we want our shape functions to be useful in p adaptive environments in 
higher dimensions. Hence, we enforce hierarchy in p, and this restricts our choice to p = 1, so the decay 
must be linear. Indeed, this is the typical and simplest choice. 



Vertex functions trace properties Linear decay Nonlinear decay 

and compatibility Hierarchy in p: ^ Hierarchy in p: ^ 


Figure 1.3: Natural vertex function compatibility and the different possible decays. 

Edge functions: In ID the edge shape functions are called edge bubbles and should vanish at the two 
endpoints of the edge (or segment). This makes them automatically compatible in ID. When p = 1 there 
are no edge bubbles, when p = 2 there is one p = 2 edge bubble, when p = 3 there is the p = 2 bubble and 
an extra p = 3 bubble giving a total of two bubbles, and so on. This results in a hierarchical construction of 
the shape functions in p. 



Figure 1.4: Potential set of of ID edge bubbles vanishing at both endpoints. 


L^. Edge functions: The ID Lf edge functions do not need to satisfy any trace properties (neither vanish¬ 
ing nor nonzero) because there is no notion of trace. They span the space of the (ID) gradients of the 
conforming shape functions, and should be hierarchical in their construction. 

1.7.2 2D 

In two dimensions, classification of shape functions is: 

: vertex, edge, and face functions, 

H (curl) : edge, and face functions, 

: face functions. 

There are two 2D elements: the quadrilateral and the triangle. Their boundaries are composed of edges and 
vertices. 














H^. Vertex functions: The vertex functions should vanish at the other (unassociated) vertices and disjoint 
edges. They should take the value 1 at the associated vertex. In 2D, at this point, this does not guarantee 
compatibility. However, dimensional hierarchy requires the (nonzero) trace of vertex functions over the 
adjacent edges to be precisely a ID vertex function (associated to the vertex in question). This results in 
full compatibility, and implies the 2D vertex functions are extensions of their ID analogues. Regarding the 
form of the shape functions across the (quadrilateral or triangle) face itself, they can have different forms 
of decay. Again, hierarchy in p will restrict our choice, so that the vertex functions we define lie in the 
lowest order space possible. Indeed, quadrilateral vertex functions present a bilinear decay, while the decay 
is linear for triangles. 



Edge functions: Edge functions should vanish at all other (unassociated) edges of the element. By 
dimensional hierarchy, for a given order p, the nonzero trace over the (associated) edge itself should take 
the form of a ID edge bubble of order p. This gives compatibility modulo edge “orientations” in 2D, 
and it implies the edge functions are extensions of the original ID edge bubbles. Now, the functions 
themselves should present a certain decay from the (associated) edge towards the rest of the element. Again, 
this choice is generally restricted by the hierarchy in p. For the quadrilateral element the decay we invoke is 
linear, but things are more complicated for the triangle element. 



Figure 1.6: Compatibility of edge functions. Depicted as well is the decay across the faces. 

Face functions: The 2D face bubbles vanish at all the edges of the element, so they are automatically 
compatible in 2D. They are constructed usually by using the edge functions previously defined (since these 
already satisfy some vanishing properties) and making some modifications to establish the remaining van¬ 
ishing conditions. Again, their construction is hierarchical in p. 

H (curl). Edge functions: The edge functions must have vanishing (tangential) trace at all other edges. 
For a given order p, the nonzero (tangential) trace over the edge itself should take the form of a ID edge 
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function of order p. This ensures compatibility modulo edge “orientations”. Now, to respect hierarchy in p, 
the edge functions themselves should present a decay (in each of their two components) which is consistent 
with the lowest order possible decay. Ultimately, the decay will be similar to that of edge functions. 

Face functions: The face bubbles have zero (tangential) trace over all edges, so they are automatically 
compatible. They are hierarchically constructed using the same ideas as for their counterparts. 

L^, Face functions: The face functions do not need to satisfy any trace properties, because there is 
no notion of trace in L^. They span the space of (2D) curls of the Tf(curl) shape functions, and their 
construction should be hierarchical in p. 

1,7.3 3D 

In three dimensions, classification of shape functions is: 

: vertex, edge, face, and interior shape functions, 

H (curl) : edge, face, and interior shape functions, 

H (div) : face, and interior shape functions, 

: interior shape functions. 

There are four 3D elements: the hexahedron, the tetrahedron, the prism and the pyramid. Their boundaries 
are composed of faces, edges, and vertices. 

H^, Vertex functions: Vertex functions have to vanish at all other (unassociated) vertices and at all disjoint 
edges and faces. Dimensional hierarchy requires the (nonzero) trace of vertex functions over the adjacent 
faces to be precisely a 2D vertex function (associated to the vertex in question). This ensures full 
compatibility of the vertex functions in 3D. 

Edge functions: Edge shape functions should vanish at all other (unassociated) edges and disjoint faces. 
By dimensional hierarchy, for a given order p, the nonzero trace over the adjacent faces should be a 2D 
edge function of order p (associated to the edge in question). Again, this ensures compatibility modulo edge 
“orientations”. 

Face functions: Face functions should vanish at all other (unassociated) faces. Dimensional hierarchy 
establishes that, for a given order, the nonzero trace over the (associated) face itself should be a 2D F[^ face 
bubble of the same order. This establishes compatibility modulo face “orientations”. 

Interior functions: The 3D interior bubbles vanish at all faces of the element, and are fully compatible 
in 3D. Our construction of these functions usually involves making some changes to the face functions 
previously defined in order to establish the remaining vanishing conditions. 
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H (curl). Edge functions: The edge functions must have vanishing (tangential) trace at all other edges 
and disjoint faces. For a given order p, the nonzero trace over the adjacent faces should be a 2D Ff(curl) 
edge function of order p. This gives compatibility modulo edge “orientations”. 

Face functions: The face functions have vanishing (tangential) trace over all other faces. Requiring the 
trace over the face itself to be a 2D Tf (curl) face bubble of a given order ensures compatibility modulo face 
“orientations”. 

Interior functions: The interior bubbles have zero (tangential) trace over all the faces and are fully 
compatible in 3D. They are constructed using the same ideas as for interior bubbles in 

H{dW), Face functions: The face functions should have vanishing (normal) trace at all other faces. The 
trace over the face itself should be a 2D face function of a given order. Compatibility modulo face 
“orientations” is then established. 

Interior functions: The interior bubbles have zero normal component at all the faces, and they are 
automatically compatible. They are constructed using the same ideas as in H^. 

L^. Interior functions: These do not have to satisfy any trace properties and they span the space of the 
(3D) divergences of the H (div) shape functions. 

1.8 Orientation Embedded Shape Functions 



Figure 1.7; Potential edge function “orientation” mismatch in 2D due to disregarding the global mesh. 

To ensure full compatibility of the shape functions along the boundaries, the concept of orientations 
needs to be introduced. The simplest example occurs in 2D, where edge functions might not match in the 
global mesh due to a simple change in coordinates over the edge itself. This “orientation” mismatch occurs 
when the edge functions in the (local) master element are transformed to the global mesh. Indeed, in the 
standard Szabo’s approach, master element shape functions are constructed with no regard for global edge 
or face coordinates, and element shape functions contributing to an edge or face basis function may not 
coincide with each other along the shared boundary. This leads to the necessity of an additional action 
during the finite element assembly process to account for local-to-global orientation changes (change of 
coordinates). Usually, these involve sign factors in the case of edges and quadrilateral faces, and more 
complicated adjustments in the case of triangle faces (see the discussion in Demkowicz et al. (2007, p.50)). 
In the context of h adaptive codes involving hanging nodes, the implementation of these modifications in 
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the assembly procedure can become quite involved, and alternative solutions to this problem are therefore 
desired. 

One such solution involves the concept of orientation embedding (Gatto and Demkowicz, 2010). Here, 
a given topological entity (an edge in 2D, or a face or edge in 3D), regardless of what elements it is adjacent 
to, is given a global orientation at the mesh level, so that it effectively owns a system of coordinates. In 
the mesh, this is equivalent to ordering the vertices in a certain order for that given topological entity. For 
instance, for an edge with vertices a and b, the vertices can be ordered as a —)■ 6 or 6 —)■ a. This choice 
defines a certain global edge orientation. This information is then passed to the master element, and the 
shape functions are defined depending on this new information. The resulting shape functions are then 
automatically compatible with each other. Naturally, at the local level, this constitutes an extension to 
the typical approach by Szabo, but at the assembly level, it simplifies fhe implemenfafion of consfrained 
approximafion (hanging nodes) by an order of magnifude. 

In view of fhese observafions, in fhis work we consfrucfed orienfafion embedded shape functions which 
lake info accounl fhe informalion regarding fhe “orienfafion” of each relevanl lopological enfily. For each 
elemenl, we explain fhese orienfafion embeddings only afler firsl presenting a complete consfrucfion of 
fhe classical (“unorienfed”) shape funclions. Hence, fhe informalion is convenienlly decoupled for ease of 
consullalion. 

1.9 Affine Coordinates 

In fhis work, we chose lo exploil simplex (barycenlric) affine coordinates lo formulaic all shape function 
conslruclions. If is well known lhaf affine coordinafes are useful when conslrucling shape funclions for fhe 
Iriangle and lelrahedron, which are simplices. However, we nole fhal wilh fhe exceplion of fhe pyramid, 
all elemenls are eilher a simplex or a Carlesian producl of simplices. Indeed, we use fhese coordinafes 
for all fhe elemenls, including fhe pyramid, where we define affine-relaled coordinafes lo complemenl fhe 
consfrucfion. 

Using affine coordinafes has many desirable advanlages. Firslly, Ihey give a solid geomelrical inluilion 
lo fhe shape funclions. Secondly, Ihey allow fhe expressions for fhe shape functions lo be used in many 
olher masler elemenl geomelries. Lasfly, Ihey play a vilal rote in fhe conlexl of orienfafion embedded 
shape funclions. Indeed, orienfafion changes are handled almosl efforllessly by simple permulafions in 
fhe argumenls of a few crucial ancillary functions (or operators). The argumenls of fhese funclions are 
precisely affine coordinafes (or affine-relaled), and Ihey are permuled in accordance lo a simple auxiliary 
permulalion function. This properly mighl be somewhal infuilive in fhe case of funclions, bul whal is 
remarkable is fhal if also holds for fhe relevanl H (curl) and H (div) funclions, where lechnically speaking, 
nonlrivial pullback maps (sometimes called Piola Iransforms) are required lo make fhese coordinale changes. 
Hence, fhese pullback maps become superfluous wilh fhe aid of ancillary operalors having affine coordinate 
funclions as Iheir argumenls. 
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Figure 1.8: Simplices in ID, 2D and 3D. 


We now define the affine eoordinates. Let vq, ... denote the vertiees of some simplex, A. Any 
point X G A ean be expressed as a eonvex eombination of the vertiees: 

N 

X = '^.SaVa. ( 1 . 10 ) 

The weights in the sum above, sq, ... ,sn, are the affine eoordinates for A. We ean think of them both as eo¬ 
ordinates in and of themselves, or funetions of the Cartesian variable x. Due to being a eonvex eombination, 
for all X G A it holds that 

N 

^Sa(x) = l, and Sa{x) > 0. ( 1 . 11 ) 

a =0 

Throughout this doeument, to ease the understanding, we shall use the following eonvention for affine 
eoordinates. 

• ID: fiQ, fii will be affine eoordinates for edges (N = 1, fia = Sa and a = 0,1). 

• 2D: Vq, Vi, V 2 will be affine eoordinates for triangles {N = 2, Va = .Sa and a = 0,1, 2). 

• 3D: Aq, Ai, A 2 , A 3 will be affine eoordinates for tetrahedra {N = 3, Xa = Sa and a = 0,1, 2 , 3). 

We will often use the following “veetor” notation for eompaetness: 

^ab (Sa; ) 'Saftc Sf,, Sc) , (L12) 

where s = ii.,v, A. Henee, for example, V 12 = ( 1 ^ 1 , ^^ 2 ) and A 031 = (Aq, A 3 , Ai). 

Explieit formulas for the affine eoordinates (in terms of Cartesian eoordinates) used for eaeh element 
will be given at the beginning of eaeh eorresponding seetion. 


1.10 Outline 

The doeument will be organized naturally starting with the simplest element in ID and then, as dimension 
and eomplexity inerease, leading into the most eomplieated elements in 3D. The order of seetions is: seg¬ 
ment (§3), quadrilateral (§4), triangle (§5), hexahedron (§6), tetrahedron (§7), prism (§8) and pyramid (§9). 
Those interested only in simplieial elements may simply read segment, triangle and tetrahedron, while those 
interested only in quadrilateral and hexadedral elements ean also skip the nonrelevant seetions. The prism 
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and pyramid sections are better appreciated after reading through all of the previous sections. As mentioned 
before, with regard to orientations, for each element there will always be a final subsection describing the 
necessary notions and modifications to implement orientation embedded shape functions. Therefore, as a 
first iteration in trying to implement the shape functions, or for those readers for which this aspect is not of 
interest, we suggest skipping those subsections. 

As a prelude to all the constructions, there is a section introducing the concept of polynomial scaling 
and the versions of Legendre and Jacobi polynomials used in our constructions. More importantly, the 
concept of homogenization is defined. This is fundamental for the elements involving triangle faces (triangle, 
tetrahedron, prism, and pyramid). 

Finally, as mentioned before, a set of tables available in Appendix E give a thorough definiton of all 
ancillary functions and shape functions presented in the text. These tables should be used as a reference by 
the reader when looking at the provided code or when implementing their own version. 

1.11 Previous Work 

Our constructions are often based either in part or in full in previous work by various collaborators in the 
field. Construction of shape functions for the quadrilateral follows Ainsworth and Coyle (2001) (see also 
Demkowicz et al. (2007)). 

For the triangle and tetrahedron, our construction is based on the concept of scaled polynomials as 
described by Karniadakis and Sherwin (1999), Schoberl and Zaglmayr (2005), Zaglmayr (2006) and the 
subsequent work of Beuchler et al. (2012a). See also Beuchler and Schoberl (2006); Beuchler and Pill- 
wein (2007); Beuchler et al. (2012b) and Beuchler et al. (2013) for more details on obtaining good sparsity 
properties via appropriate selection of Jacobi polynomials. Contrary to their work, here we study the classi¬ 
cal H (curl) and H (div) conforming Nedelec and Raviart-Thomas spaces having the property that they are 
affine invariant, and being compatible (at the space level) with the spaces proposed for the pyramid. Other 
interesting shape functions for the tetrahedron include those of Ainsworth et al. (2011) based on Bernstein 
polynomials. Lastly, it is worth noting that Zaglmayr (2006) also presents a unified construction of the 
hexahedron and prism to complement the tetrahedron, but does not include the pyramid. 

The prism element is a Cartesian product of the 2D triangle and ID segment. The prism shape functions 
therefore utilize constructs from the triangle and segment. 

Construction of pyramid shape functions builds on the fundamental work of Nigam and Phillips (2012) 
and their first family of pyramid spaces. The spaces are natural for (parallelogram-based) affine pyramids, 
but as evidenced by Bergot and Durufle (2013), they also have other attractive properties in a non-affine 
setting. We also note that Bergot et al. (2010) and Bergot and Durufle (2013) have contributed to the work 
on higher order pyramid shape functions, but their spaces and shape functions are different. 

The idea of orientation embedded shape functions follows the work of Gatto and Demkowicz (2010) 
and stems from discussions with Joachim Schoberl dating back to the Vienna WCCM congress in 2002. 
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2 Polynomials Prelude 


2.1 Notation 

The polynomials of order p with arguments x G M will be denoted by 

V^{x) = spanjx-^ : j = 0,... ,p} . ( 2 . 1 ) 

Similarly, in two dimensions the polynomials of total order p with arguments (x, y) G are denoted by 

V^{x, y) = span{x*y '■ i > 0, j > 0,n = i + j < p} , (2.2) 

while the homogeneous polynomials of total order p are denoted by 

V^{x, y) = spanlx^y-^ : i > 0, j > 0, i + j = p} . (2.3) 

Similar definitions apply to polynomials of three variables. Moreover, when the domain is elear from the 
eontext, we will simply refer to V^{x, y) and V^{x, y) as and respectively. 

Define 

Q^’'^(x, y) = 'P^(x) (g) 'P^(y) = span{x*y^ : 0 < i < p, 0 < j < g} , (2.4) 

and similarly for Q^’'^’'’(x, y, z). When fhe variables are clear from fhe confexf, fhese spaces are simply 

wriffen as and QP’'?’'" respectively. 

The nofafion for vecfor valued polynomial spaces will be 

(pp )2 = pp X pP , ( 2 . 5 ) 

and similarly for and fhe vector valued homogeneous polynomials, {V^)^, N = 2,3. 

2.2 Scaled Polynomials 

Given an order i univariate polynomial, 'pi{x) G 'P*(x), we define fhe corresponding scaled polynomial 

V'i(x;f) = • (2.6) 

Obviously, ipi{x] 1) = p’i{x), so fhe scaled polynomials define fwo variable polynomial extensions into fhe 
(x, t) space. Furfhermore, fhe reader may observe fhaf 'ipi{x; t) is homogenous of order i as a polynomial in 
fhis space, i.e. V'i(x; t) G T’*(x, t). 

2.3 Legendre Polynomials 

In fhis work, we will use Legendre and Jacobi polynomials for fhe consfrucfion of all shape functions. 
Should fhe reader wish fo work wifh differenl families of polynomials, our consfrucfion easily generalizes 
as discussed in Appendix A. 
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The classical Legendre polynomials comprise a specific orthogonal basis for L^(—1,1). Truncated to 
the first p elements, {Pi : i = 0,... ,p — 1}, the Legendre polynomials^ are a basis for the space of (single 
variable) polynomials of order p — 1. 

Of many properties of the Legendre polynomials, we list the following recursion formula 

h{y) = 1, 

Pi{y)=y, (2.7) 

iPiiy) = (2i - 'i-)yPi-i{y) - {i- i)Pi- 2 {y), for t > 2 , 
and the derivative formula for i > 1, 

(2z + l)Piiy) = ^ (^Pi+i{y) - Pi-iiy)^ , (2.8) 

which are both well known in the literature. We also make note of the orthogonality relationship 

j Pi{y)Pj{y)(iy = 0, (2.9) 


This relationship, and the definition Po{y) = 1, leads to the zero average property 


Pi{y)dy = 0, for t > 1. 


'-r 


( 2 . 10 ) 


Shifting. The range of affine coordinates is always [0,1] (see (1.11)). Although not clear at the moment, 
this implies that we want to have the zero average property over the interval [0,1] instead of [—1,1]. We can 
obtain this property by composing each Legendre polynomial above with the shifting operation 

y^2x-l. ( 2 . 11 ) 

The (shifted) Legendre polynomials over [0,1] are defined for t > 0, 

Pi{x) = Pi i2x - 1) . (2.12) 


Scaling. The (shifted) scaled Legendre polynomials are defined by (2.6) as 

Piix; t) = P, (I) f = P, (2 (I) - 1 ) f = Pi{2x (2.13) 

where the domain in x is [ 0 , t]. 

"'Although the common notation for the classical Legendre polynomials is Pi, we choose to denote the elements in this set with 
~ as we will only need this definition temporarily. 
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This set of scaled Legendre polynomials obeys a recursion formula similar to (2.7): 


Po{x;t) = 1, 

Pi{x; t) = 2x — t, (2.14) 

iPi{x; t) = {2i — l)(2x — f)Pj_i(x; t) — {i — l)f^Pj_ 2 (x; t), for i > 2 . 


Moreover, we carry over the orthogonality and zero average properties from (2.9) and (2.10) to the scaled 
domain [0 ,f]: 

/ Pi{x]t)Pj{x;t) dx = 0 , if i / j , 

Jo 

/ Pi{x; t) dx = 0 , fori>l. 

Jo 


(2.15) 


Integrated Legendre Polynomials. For all i > 1, we define the (scaled) integrated Legendre polynomials, 

px 

Li{x;t) = / Pi-i{P,t)dx, (2.16) 

Jo 

where of course Li{x) = Li{x; 1). Notice that pP{x) = span({l} U {Li : f = 1,. .. ,p}). By construction, 
the Li are seen as elements of and as a result, their pointwise evaluation is understood to be well defined. 
Therefore, recalling fhe zero average properly of Ihe Legendre polynomials, we observe fhal, 

Lii0) = Li{0-t) = 0 = Li{f,t) = Liil), fori >2. (2.17) 


Nexl, (2.8) mofivales fhe formulas for computing fhe inlegraled Legendre polynomials: 

Li(x; t) = x , 

2{2i — l)Li{x; t) = Pi{x; t) — t), for i > 2 . 


(2.18) 


Clearly, 

^Li{x;t) = Pi-i{x;t). (2.19) 

ox 

Derivatives of Li(x; t) wifh respecf lo t will also be necessary in our compulations. For Ihis, we define 

Ri{x) = {i + l)Lj+i(x) — xPi{x), for i > 0 , (2.20) 

which fhe reader may observe is an order i polynomial. In Appendix A we show fhal 

^Li{x]t) = Ri-i{x;t). (2.21) 

Obviously Ro{x; t) = 0, and by use of (2.14) and (2.18), one can reduce (2.20) lo 

Ri{x-,t) =-^(^Pi{x-,t)+tPi-i{x-,t)'j , forf>l. (2.22) 
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2.4 Jacobi Polynomials 


Motivated by Beuchler and Schdberl (2006) and Beuchler and Pillwein (2007), we use Jaeobi polynomials in 
our eonstruetions of elements involving triangle faees. The (shifted to [0,1]) Jaeobi polynomials, P- ' 
a,f3 > —1, form a two parameter family of polynomials ineluding the Legendre polynomials previously 
defined = p^)_ Jaeobi polynomials have similar reeursion formulas as the Legendre polynomials. 

One may find a seleefion of sueh formulas in Beuehler and Pillwein (2007). For our purposes, we will only 
eonsider fhe ease /3 = 0, so fhaf from now on 

Jaeobi polynomials are also orfhogonal in a weighfed spaee. Assuming fhe sealing operation dis- 
eussed previously, we have fhe orfhogonalify relafion 

x“P"(x; t)P^{x; f) dx = 0 , ifi ^ j , (2.23) 

whieh for a / 0 no longer implies fhe zero average properly. 

The following is fhe reeursion formula we use lo eompule fhe [0, t] Jaeobi polynomials: 

P^{x-t) = l, 

Pi{x;t) = 2x — t + ax, (2.24) 

aiP^{x; t) = bi [ci{2x — t) + a^t) P^_i{x; t) — dit^P^_ 2 {x] t), for i > 2 , 

where 

tti = 2i{i P a) (2i + a — 2), 
bi = 2i P a — 1, 

Ci — {2i T cr) (2i P cx — 2) , 
di = 2{i P a — 1) (f — 1) (2i + a) . 

We remark lhal ofher reeursive relations in weighl and order lo eompule Jaeobi polynomials, sueh as 
(a + z)P“(x; t) = (a + 2i)P““^(x; t)pitP^_i{x‘, t), were experimenlally found lo be numerieally unslable 
as eompared lo fixing a value of a and using (2.24), so lhal Ihe laller approaeh is reeommended. 

Integrated Jacobi Polynomials. Finally, we define the (sealed) integrated Jaeobi polynomials for i > 1: 

Lf{x;t)=f 7^“i(x;f)dx, (2.25) 

Jo 

with Lf{x) = Lf{x; 1). Note that beeause of the absenee of the zero average property, we eannot deduee 
that Lf{l) = 0, and in general, this does not hold. However, it is obvious that for all a > —1, 

L“(0) = Lf(0; t) = 0, for i > 1. (2.26) 
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We evaluate the integrated Jaeobi polynomials using the following relations:^ 

t) = X , 

Lf{x; t) = aiP^{x, t) + bitP^i(x; t) - Cit^P^_ 2 {x] t) , for i > 2 , 

where 


^ + a 


a,' = 


h = 


{2i + a — l)(2i + a) 
a 

{2i + a — 2){2i + a) 
i - 1 


Ci =- 


{2i + a — 2)(2i + a — 1) 
As in the case of the integrated Legendre polynomials, we find that 


■^Lf{x-,t) = ^Lf{x;t) = Rf_^{x-,t), 


where again, 


Rf{x) = {i + l)L“+i(x) - xP°‘{x). 

Obviously Rq{x; t) = 0, and by use of (2.24) and (2.27), one can reduce (2.29) to^ 

Kix-R) = - (pt{x]t) + tP^_i{x]t)) , for i > 1 

2^ + a V / 

2.5 Homogenization 

Definition. For an order i polynomial 

Tpi G V^Sl, . . . , Sd) , 

we define the operation of homogenization [of order i) as a linear transformation 

[ • ] : V\si, ...,Sd) —^ Vfso, si,..., Sd), 

where 


[V’i](so,si, ...,Sd) = fi 


Si 


Sd 


So + • • • + Srf ’ ’ So + • • • + Sd 


(so + • • • + Sdf 


(2.27) 


(2.28) 

(2.29) 

(2.30) 


(2.31) 


Notice that homogenization is a form of scaling, and as such, it is forming an extension of the particular 
case in which sq + si • • • + Sd = 1. It is not a coincidence that this is precisely the property that affine 
coordinafes safisfy (see (1.11) in §1.9). Moreover, note fhaf [fiii] is always a homogeneous polynomial of 
degree i, so we have fhe following scaling properly for all scalars 7, 


[V’i](7So,7Si,- • • ,7Sd) = 7 *[V’i](so,si, •.., Sd) • 


(2.32) 


^cf. (2.9) in Beuchler and Pillwein (2007). 
®cf. (2.16) in Beuchler and Pillwein (2007). 
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One will observe that for the partieular ease d = 1, 

[V’iKsoi'Si) = V'i('Si;so + Si) • (2.33) 

Therefore, we see that 

[tpi]{so,si) = ipiisi;!) = ^pi{sl) , if So + ■51 = 1, (2.34) 

where we remind the reader that the ID affine eoordinates satisfy preeisely this property. 

Moreover, take the case of the integrated Legendre polynomials and recall property (2.17). It follows 
that for all So, Si, 


[Li](so,0) = Li(0;so) = 0 = Li(si;si) = [Li](0,si), for i > 2 . (2.35) 


In the d = 2 case, one can make another useful observation. Let Xj be a one variable polynomial of 
order j, and consider the homogenization of the i + j order (two variable) polynomial 


= [^pi]{so, Si ) Xj{l - So - Si ) . 


(2.36) 


In this case, using (2.32), we find 


[V’ij](.S2,So,Si) = [V'j] 


So 


Si 


So + Si + S 2 So + Si + S 2 

■j['tpi]{so,si)xj 


)v(- 


So + Si 


(so + Si + S 2 )*’^-^ 


(so + Si + S 2 ) 

= [V'i](so,si)[xi](so + Sl , S2 ) . 


S2 


So + Si + S2 


So + Si + S 2 . 

(so + Si + 52 )*^-^ (2.37) 


This inspires fhe following definition for single variable polynomials, V’t and Xj, of order i and j respec¬ 
tively: 

[V'i,Xi](so,si,S 2 ) = [V’i](so,si)[xj](so + S1 , S2 ), ( 2 . 38 ) 

where if is clear [ipi, Xj] £ 'P*"'“-^ (so, si, S 2 ) is a homogeneous polynomial of order i + j. 

Again, observe fhaf 


[V'LXi](so,si,S 2 ) = [V’i](so,si)xi(s 2 ), if so + Si + S 2 = 1, (2.39) 


where we remind fhe reader fhaf fhe 2D affine coordinafes salisfy precisely fhis properly. 

Finally, using properlies (2.17) and (2.26) of fhe inlegraled Legendre and Jacobi polynomials, if follows 
fhaf for all so, si, S 2 , 


[Li,L“](so,si,0) = [Li,L"](so, 0 ,S 2 ) 


[Li,L“]( 0 ,si,S 2 ) = 0, forf>2,j>l. 


(2.40) 
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3 Segment 


0 * *1 

Figure 3.1: Master segment with numbered vertices. 

The ID simplex is the segment or edge. The master element is defined as the unit interval (0,1), and it 
is illustrated in Figure 3.1 with a parameterization given by ee[0, 1]. 

Denote vertex a by Vq. The definition of affine eoordinates (see (1.10)) sfafeslhaf^ = ^o(0^o+Fi(0^i’ 
wifh //o(0 + Fi(C) = Fo(0 ^ 0 ;Ui(^) > 0 for all ^ G [0,1]. Henee, /jq is fhe weighl related fo vq 

and Hi is fhe weigh! relafed fo t;i. For our masfer elemenf, wq = 0 and t;i = 1. If fhen follows fhaf fhe ID 
affine eoordinates for fhe segmenf, ho, Ft’ the mosf basie linear funelions. They are wriffen explieifly 
below for our masfer elemenf: 


Fo(e) = l-C, Fi(6 = e- (3-1) 

The gradienfs of fhe affine eoordinafes (in ID) are 

V;uo(0 = -1, Vmi(0 = 1- (3-2) 


Exact Sequence 

The ID exaef sequenee ean be eonsulfed in §1.2. Polynomial spaees are subsefs of bofh and (in 
(0,1)), so we may eonsider a fruneafed polynomial approximafion (of order p) fo and induee a new 
(diserefe) polynomial exaef sequenee 

pp pp-i ^ (3.3) 

where = 'Pt’(^). in the notation of (1.8), and = V^~^. 

3.1 Shape Functions 

The set of all shape funetions defined in fhis seefion will form a basis for fhe spaee whieh has dimension 
p + 1. In fael, fhere will be 2 verfex shape funelions and p — 1 edge shape funelions. They will all be linearly 
independenf and be eonlained in so Ihey will elearly form fhe desired basis. 

3.1.1 Vertices 

As previously mentioned, eaeh vertex is linked to an affine eoordinate. For instanee, vq is linked to hq. It is 
then quite natural to eonsider the affine eoordinate itself as the associated vertex shape funetion to vq\ 

4>^{0 = doiO ■ 
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Indeed, it satisfies all the desired traee properties, sinee it takes the value 1 at no = 0, and 0 at ni = 1. 
Moreover, it deeays linearly to the other vertex, so that it lies mV^, and respeets the hierarehy. Having 
the vertex funetion of the form instead, would give a (faster) nonlinear deeay, but then the funetion 

would be dependent on p and the hierarehy would be broken. 

In general, the vertex funetions, along with their gradients are, 

= V/ra(0 , (3.4) 

for a = 0,1. Clearly, there are a total of 2 vertex funetions (one assoeiated to eaeh vertex). 

3.1.2 Edge Bubbles 

Reeall from §2.3 that the Legendre polynomials are seen as elements of whieh have the zero average 
property (over [0,1]), so that the integrated Legendre polynomials (of order 2 and higher) are elements of 
whieh vanish at 0 and 1 (see (2.17)). These are preeisely the desired eharaeteristies for edge bubbles. 
Henee, the edge funetions are defined as: 

YiiO = i = 2,...,p. 

This formula is perfeetly valid and quite simple. However, along this doeument, the use of affine eoordinates 
will be enforeed as mueh as possible. The reasons for this will beeome elear as we move into higher 
dimensions. Indeed, notiee that due to pi{C) = one ean write Lj(C) = Moreover, sinee 

IJ-o + = 1. by (2.34) it follows 


Li(/ii) = Li{pi] 1) = [Lilipo, pi ). 


With this in mind, eonsider the following more general setting. 

Definition. Let sq si be arbitrary functions of some spatial variable in with N = 1,2, 3. Denote 
by ps the order in the coordinate pair (sq, si)- Then 

())f (so, Si) = [Li]{so, Si) = Lj(si; So + si), (3.5) 


for i = 2,... ,ps. The gradients, understood in are 


V(/)f(so,Si) = [Pi_i](so,Si)Vsi + [i?j_i](so,Si)V(si + Sq) 

= Pi_i(si; So + si)Vsi + Ri-i{si; so + si)V(si + so). 


(3.6) 


Clearly, the definition of (f^, involving homogenization, ean be thought of as an extension of our more 
simple ease. This is the first of the so-ealled ancillary functions whieh are defined in this work. It is 
highlighted as an important definition, beeause it will be used multiple times throughout the text in more 
general settings. Here, the superseript E stands for edge, and one should think of this topologieal entity 
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when looking at this function. Also, note its arguments, sq, si, are meant to be affine coordinates (or at least 
affine-related). 

Rewriting (2.35), it follows that for any sq, si, and all i >2, 

4,f{0,si) = 4>f{so,0) = 0- (3.7) 

As observed, when the coordinates are ID affine coordinafes, like in fhis case, fhe formulas for (j)f and 
ifs gradienf are simplified. For fhis, record fhe nexf remark. 

Remark. Lei = 1 — ni, where /xi is an arbitrary function of some spatial variable in N = 1,2, 3, 

and where p is the order in the coordinates {po, pi). Then for alii = 2,... ,p, 

4'f{lto,Ti) = Li{pi), V(l)f{po,pi) = Pi-i{pi)Vpi. (3.8) 

From now on, shape funcfions will be wriffen in terms of ancillary funclions and fhe affine coordinafes 
of fhe elemenf being analyzed. Indeed, all fhaf is required fo compute ff and Vff are sQ) •si and Vsq, Vsi 
(since we already know from §2.3 how fo compufe fhe scaled versions of Li, Pi and Ri). For fhe segmenf, 
■So = Po, Si = pi, so fhis information is in (3.1) and (3.2). 

Af firsl, fhis approach mighf seem fo be overcomplicafed given fhe simplicify of fhe inifial formula 
(which does nol involve scaling). However, compufafionally speaking, fhis mofivafes coding ff and Vff, 
which will be observed fo be fundamenlal as fhe documenf progresses. If desired, wifhin fhe ff subroufine, 
one could decide fo separafely handle fhe special sifuafion where sq, si are ID affine coordinafes, in which 
case fhe simplificafion shown in (3.8) would fhen hold. More imporfanfly, when orienfafions become rele- 
vanf, fhey will be handled fhrough permufafions of fhe argumenfs of ff. Hence, having everyfhing wriffen 
in ferms of ff, and in general, in terms of ancillary funclions, is highly desirable. 

Recalling fhaf /Tot = (fO) Ft)’ the shape functions and fheir gradienfs are fhen defined as 

m) = </>f (foi(6) , ^m) = v</.f (Moi(e)), ( 3 . 9 ) 

for i = 2 ,... ,p. There are p — 1 edge bubbles for fhe segmenf. 

As previously mentioned, if is clear fhaf 0|(O) = <?i>|(l) = 0, by fhe properfies of fhe infegrafed Legendre 
polynomials when i >2 (see (2.35) or (3.7)). Hence fhe frace properties are satisfied. 

3.2 Shape Functions 

The colleclion of conforming shape funclions is simple and motivated from fhe exacl sequence. In ID, 
all funclions are realized as gradienfs of functions. To resemble Ibis properly al fhe discrete level, we 
simply consider fhe linearly independenl derivalives of f'' and 0?. Clearly fhey will be a basis for P^~^. 
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3.2.1 ^2 Ejggg 

The edge shape funetions are the Legendre polynomials, which written in terms of affine coordinates are 

, o.io) 

for i = 0,... ,p — 1. There are p such edge functions and they span V^~^. The apparently trivial fac¬ 
tor V Pi (^) makes the expression coordinate free, so it takes the same form independent of any (possibly 
nonlinear) transformations. 

3.3 Orientations 

In ID, the trace is simply the two endpoints of each element, and it is clear that shape functions of adjacent 
elements will have full compatibility at the vertices. However, in 2D and 3D, the boundaries involve edges 
and faces. Achieving this compatibility is nontrivial. By dimensional hierarchy, edge functions of elements 
in higher dimensions will involve (through the trace operation) the ID edge functions defined in fhis section 
(see §1.5). In view of fhis, if is nafural fo explain edge orienfafions af fhis time. 

3.3.1 Edge Orientations Explained 



In local edge coordinate when o = 0 In local edge coordinate when o = 1 


Figure 3.2: Edge orientations. 

In 2D, if one naively disregards how the elements are placed in the global mesh, and proceeds to define 
all shape funclions af fhe (local) master elemenf level, one mighf end up wifh shape funclions fhaf, when 
Iransformed back info fhe original mesh, are incompafible across shared edges (see Figure 1.7). Wifh ori- 
enfafion embedded shape funclions fhis problem is avoided by faking info accounl more informalion of fhe 
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global mesh. This is done by giving each mesh edge its own global orientation, and is represented by a 
global coordinate or equivalently by a global edge vertex-ordering. For example, given an edge at the 
mesh with vertices a and b, a global edge vertex-ordering of the form a —)■ 6 means that has its origin at 
a and points from a to b. This information is then passed to the particular master element, where the edge 
has its own fixed local orientation, represented by the local coordinate or equivalently by the the fixed 
local edge verfex-ordering of fhe form (nofe fhe dashed arrow for local orderings). Viewed af fhe 

local level, fhe global coordinafe can eifher coincide wifh fhe local coordinate ^ or poinf in fhe opposite 
direcfion. To reflecf fhese two possibilifies, fhe orienfafion parameter o is infroduced. If o = 0, fhis means 
fhe local and global coordinafes coincide, and ofherwise o = 1. All fhis is depicfed in Figure 3.2. 

Due fo fhe use of affine coordinafes and fhe form of fhe ancillary funclions proposed in fhis work, fhese 
orienfafion problems can be readily fackled. To ensure full compafibilify, we wanf fhe shape funclions over 
a given edge lo be immovable when observed in fhe global coordinafes (like in Figure 3.2). This is achieved 
by evaluating fhe ancillary functions wifh global coordinafes. Unforlunafely, fhe available coordinafes pro¬ 
duced by fhe master elemenl are fhe (fixed) local coordinafes. Hence, fhe idea is lo apply a local-to-global 
transformation over fhe edge, which will obviously depend on fhe orienfafion parameler o. Such a Iransfor- 
mafion is complelely nalural in fhe conlexl of affine coordinates, since fhis only involves permulafions of 
fhese coordinafes. Indeed, a simple permulalion function dependenl on o, denoted by a^, will represenl fhis 
Iransformalion. 


Definition. Let sq and si be arbitrary variables, and let o = 0,1 be the edge orientation parameter. The 
edge orientation permutation function, is defined as 


(sO,'Sl) 


cr^(so,si) = (so, Si) 
Crf('S0,Sl) = (si,So) 


if 0 = 0, 
if 0 = 1. 


(3.11) 


To explain fhe definition of a^, note lhal in (3.11), if one links sq to the local verlex [o] and si lo fhe 
local verlex QJ, Ihen fhe locally ordered pair (sq, si) represenls fhe local coordinafes. If is ordered in fhe 
sense lhal so comes firsl and si comes second, and fhis is meanl lo correspond wifh fhe fixed local ordering 
[0]--)-|l], where comes firsl and 0 comes second. Similarly, Ihere are globally ordered pairs which 
depend on fhe parameter o. Indeed, in Figure 3.2, looking al fhe global edge vertex-ordering a^b, Ihere is 
an induced global verfex-ordering of fhe Iwo verlices [O] and Q]. If is [()]—)• Q] if o = 0, and )■ [o] if o = 1. 
Hence, fhe global coordinates are represented by fhe globally ordered pairs (so, si) if o = 0 and (si, so) if 
0 =1. Therefore, in fhis sense, af is a local-lo-global Iransformalion. 

Now, all lhal is required is lo compose fhe edge ancillary funclions and Iheir differential form (Ihose 
wifh superscripl E) wifh fhis permulalion funclion af. Thus, in 2D and 3D, all inslances of ff and Vff in 
fhe shape funclions should be replaced wifh ff o af and Vff oaf respeclively. The resulling funclions are 
Ihen said lo be orientation embedded shape funclions. More concrele examples will be given in fhe 2D and 
3D elemenls as fhe documenl progresses. 


25 



4 Quadrilateral 


?2 


4 


3 


1 2 

Figure 4.1: Master quadrilateral with numbered vertices. 

The master element for quadrilaterals, which is (0,1)^, is shown in Figure 4.1 in ^ = (^ 1 ,^ 2 ) space. 
The master quadrilateral is clearly a Cartesian product of two segments. 

Due to the product structure, there are two pairs of ID affine coordinates: 

Fo(6) = 1-?i, Fi( 6)=6 =► V/ro(6) = ("qM , = 

^ ^ ^ ^ (4.1) 

/^o( 6 ) = i- 6 ) fi(6) = 6 =► v^o( 6 ) = ) v^i(^ 2 ) = (^?) • 

These will be used explicitly or implicitly in the formulas that follow. 

Again, vertex a is denoted by Va, so that vi = (0,0), V 2 = (1, 0), V 3 = (1,1) and = (0,1). These 
vertices are related to the affine coordinafes jusf as in ID. For example, over edge 12 (or edge 43), ^o(?i) 
is fhe weighf relafed fo vi, while is the weight related to V 2 (and similarly with V 4 and ^ 3 ). Indeed, 

given a point (Ci,0) on edge 12, it holds that (^i,0) = + fJ-i{^i)v 2 - The same way, /io(? 2 ) is 

related to vi in edge 14, so that vi is linked to both /xo(Ci) and /Uo(^ 2 )- A similar assertion holds for each 
vertex, where each is linked to two affine coordinafes. Now, looking af fhe edges, nofe fhaf ^ 10 (^ 2 ) takes the 
value 1 over edge 12 and 0 at opposite edge 43. This way, each edge is linked to one affine coordinafe. 


Exact Sequence 

Recall fhe 2D exacf sequence for simply connecfed domains (1.2) and ifs rofafed analogue (1.5). The corre¬ 
sponding polynomial exacf sequences are 

QP,q gp-hg ^ gP-'i'-i ZA^gp-1,9-1 

’ (4.2) 

gp ,9 TlTi^gp, 9 -i x gP-h 9 ^gp-1,9-1 ^ 

where = Q^’'^{Ci,C 2 ) = 'P^iCi) ® ’P'^(C 2 )- These are the standard Nedelec’s spaces (1980) of the 
first type for the quadrilateral. Note here the natural anisotropy of the element, which has order p in the 
direction and a potentially different q in the ^2 direction. The hierarchy should be maintained in both p and 
q separately. This is associated to the notion of local p adaptivity. It will sometimes be convenient to refer 
to Pa as the order in the direction, so that pi = p and p 2 = q. 
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4.1 Shape Functions 

It will be clear that all shape functions lie in and that they span the space. For this, one will only require 
the linear independence of the shape functions, which will be evident, and a judicious count of them, which 
will give {p + l)(q + 1) (the dimension of 

Also, note that due to the Cartesian product structure, there is a natural separation of variables, and one 
expects the shape functions to be tensor products of the relevant ID functions for the edges and vertices. 
That is, tensor products of Pai^b) and for a = 0,1 and 6=1,2. Fortunately, this is the case. 

4.1.1 Vertices 

As mentioned before, each vertex is linked with two affine coordinates, and the associated vertex function 
is precisely the tensor product of these two coordinates. For instance, vi is linked to po{Ci) and ^ 0 (^ 2 ), so 
its associated vertex function is 

= Fo(6)fo(6) • 

It satisfies all fhe desired properties, since if vanishes af fhe disjoinl edges 23 and 34, and more imporfanfly, 
ifs frace over fhe adjacenf edges is a ID verfex function associated fo fhe vertex. For insfance, over edge 
12, where po{C 2 ) = 1. its frace is /ro(?i)? which is fhe ID verfex function associated fo vi over fhe edge 
12. Finally, fhe function decays bilinearly and is in fhe lowesf order possible space, so fhaf if respecfs 
fhe hierarchy in bofh p and q. 

More generally, fhe verfex funcfions and fheir gradienfs are, 

= /^a(6)Ffe(6), V</)"(0 = /ra(6)V//;,(6) + F6(6)VMa(6) • (4.3) 

for a = 0,1 and 6 = 0,1. There is a fofal of 4 verfex funcfions (one for each verfex). 

4.1.2 Edges 

To ease fhe comprehension, fake for example edge 12, where .^2 = 0. The idea for fhe edge funcfions is 
fo use fhe segmenf bubbles (ff in ^1 and blend fhem wifh a linear funclion in ^ 2 - Thai is, blend fhem wifh 
fhe linked ID affine coordinafe, /Uo(^ 2 )- Hence, fhe associaled edge funcfions will be fhe fensor producls of 

(Pfimiii)) and poii-2), 

(PiiO = Fo(6)</>f (foi( 6)) = (1 - 6)Ti(6), 

wifh i = 2, ... ,p. The frace properties are salisfied mainly due fo fhe vanishing conditions of fhe Lj af 
fhe endpoinls (see (2.17)), which are resfafed in terms of (pf in (3.7). The vanishing properties are easily 
observed in fhe simplified form, (1 — ^ 2 )Ti(^i), bul we will write fhese Iraces in terms of the ancillary 
functions and the affine coordinafes. For Ihis, if is useful fo wrile fhe boundary reslricfions in lerms of affine 


27 



coordinates. For example, over edge 12, which has equation ^2 = 0, /2oi(?2) = (^ 0 (^ 2 ), ^ 1 (^ 2 )) = (Ij 0). 
Using these natural relationships, the edge traces are 

</’l(6l6=o = Fo(6)<?^*^(Foi(6))l/toi(6)=(i,o) = 1 • (pfiP-oiiCi)) = , 

</’l(0ki=i = Fo(6)<?^’f (Foi(?i))|/tQi(gi)=(o,i) = Fo(6)</>f (0,1) = 0, 

</’l(0k2=i = Fo(6)<?^’f(Foi(?i))|/toi(€2)={o,i) = O' 4>f{Poi{Ci)) = 0, 

(/>|(Olci=o = ^o(6)?^*^(Foi(6))l/toi(€i)=(i,o) = Fo(6)'/>f (1,0) = 0. 

Hence, the desired vanishing properties are satisfied, and more importantly, over the edge 12 itself, the trace 
is (/)f (/loi(Ci)) which as expected is a ID edge bubble. Note here the decay towards the rest of the 
element, represented by the blending function ^ 0 (^ 2 ), is linear. Indeed, the edge functions for edge 12 lie in 
and they respect the hierarchy in q. 

Next, we will give a geometrical representation of the edge 12 shape functions presented above. Recall 
that over edge 12, vi is linked to /io(?i), while V 2 is linked to Hence, /loi(^i) = (fio(?i), Fi('^i)) 

actually represents a point in the edge, 

/^o(6)(o) +Fi(6)(o) = ( 0 ) ■ 

This can be interpreted as a projection to edge 12 from an arbitrary point (^ 1 , ^ 2 ), 


(6,6) 


( 6 , 0 ). 


The trivial projection consists simply of finding fhe infersecfion P' = (6,0) of fhe edge wifh fhe normal 
projecfing line passing fhrough fhe original poinf P = (6,6)- H is heller illuslraled in Figure 4.2. 



blend 



Figure 4.2: Edge projection from P to P', and the logic project —evaluate —> blend. 

After the original point is projected to the desired edge, it is evaluated at that edge, and finally if is 


blended linearly: 


'/'KO = Fo ij. 2 ) 4>i ( fo66) ) = (1 - 6) 6(^6^) • 

blend project blend project 


evaluate evaluate 

This blending represenls and exlension or lifling lo fhe resl of fhe elemenl. The whole process of 


projecfing 


evalualing 


blending 


28 













is extremely important, sinee it is used in the construetion of shape funetions for all remaining elements. 
Moreover, as depicted in Figure 4.2, it gives a geometrical interpretation to the formulas. It should be 
mentioned that if orientations are to be handled, they are taken care of only at the level of evaluating, where 
a local-to-global transformation will need to be prepended to the original evaluating procedure. Projecting 
and blending are unaffected by orientations. 

The general formula for edge functions is 

= Fc(4)</>f (FOl(Ca)) , vm) = Fc(4)V</.f (/2oi(ea)) + , (4.4) 

where i = 2,... ,pa, (a, b) = (1, 2), (2,1) and c = 0,1, with pa being the order in the coordinate. For 
example for edge 12 (linked to pdCb) = this would correspond to (a, 6) = (1,2), c = 0 and 

Pa = P, for edge 23 (linked to /ri(^i)) it is (a, b) = (2,1), c = 1 andp^ = q, and so on. For each edge there 
are Pa — 1 shape functions, leading to a total of 2 {p — 1) + 2 {q — 1) edge shape functions. 

4,1,3 Face Bubbles 

The quadrilateral face bubbles can be naturally defined as the tensor product of ID edge bubbles, 

for i = 2,... ,p and j = 2,... ,q. Using (3.7), it is clear the vanishing conditions over all four edges are 
satisfied. This motivates a more general definition. 

Definition, Let (sq, si) and (to, ti) be two pairs of coordinates which are arbitrary functions of some spatial 
variable in N = 2,3. Let ps be the order in the (sq, si) coordinates, andpt be the order in the {to, ti) 
coordinates. Then 

= </>f('So,si)0j (fo,fi), (4.5) 

for i = 2,... ,ps and j = 2,... ,pt. The gradients, understood in are 

^4’?j{so,si,to,ti) = (l)f{.so,si)V(j)f{to,ti) + (/)^(fo,fi)V(/)f(so,si) • (4.6) 

Rewritten in terms of fPj, the general formulas for the bubbles and their gradients are, 

= 4>?jimi^i),mi^2)), v4(0 = v</.°(/7oi(6),m(6)), (4.7) 

where i = 2,... ,p and j = 2,... ,q. There are a total of (p — 1) (q — 1) such functions. 

4,2 i?(curl) Shape Functions 

It will be clear that all shape functions will lie in x Moreover, after all functions are defined, 

a rigorous count will give p{q + 1) + (p + l)q, which is precisely the dimension of QP“^’'? x QP’'^“^, so that 
the shape functions span the desired space. 
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4.2,1 H (curl) Edges 


First, take for instance edge 12. As mentioned in §1.6, the tangential trace of the Ff(curl) edge func¬ 
tions should be a ID shape function. From (3.10), the ID edge functions (with coordinate .^i) are 
[-Pi](Foi(^i))- Meanwhile, note the tangential vector to edge 12 is (1,0) = V When coupled with a 
blending factor, ^ 0 (^ 2 )> representing a linear decay (like that of H^), this suggests, 

Em = Mo(6)[^’*](m(6))v^i(a) = (1 -6)^*(a)(j), 

for i = 0,... ,p — 1. Next, the trace properties are checked. For this, note that (0,1) = V/ii(^ 2 ) is the 
tangent direction to the edges 23 and 14 (where = 1 and = 0 respectively). Hence, 

tr(Ef(e))ie=o = Emmi,)=(m ■ iv2 -v,) = i- mmim • i = [^*](Foi(ei)), 

tr(E'r( 6 )ki=i = ^i'(Olm(«i)=( 0 ,i) • (^3 - V 2 ) = /ro(6)[-Pi](0,1) • 0 = 0, 

tr(E'r( 6 )l 6 =i = ^i'( 0 lm( 6 )=(o,i) • (^3 - ^ 4 ) = 0 • [Fi](/roi( 6 )) -1 = 0 , 

= ^i'(OI/Joi(a)=(i,o) • (^4 - ^^i) = /ro( 6 )[-Pi](l> 0 ) -0 = 0 . 

The trace properties are then satisfied. Inspired by first order Whitney functions, this motivates the following 
more general definition. 

Definition. Let sq si be arbitrary functions of some spatial variable in with N = 2,3. Denote by 
Ps the order in the coordinate pair (sq, si). Then 

Ef (so, si) = [Pi\{so, si)(soVsi - siVso), (4.8) 

for i = 0,... ,ps — 1, and where the gradients are understood in M^. The curls are 

V X Ef(so,si) = (f+ 2)[Pj](so,si)Vso x Vsi. (4.9) 

Here, if = 2, the curl and cross product take the form described in (1.3). Note Ef involves the 
gradients of its entries, so that it is actually a differential operator assumed to be acting on a functional space 
(the entries sq, si functions). Hence, the use of the term ancillary operator is perhaps more appropriate in 
this case. The final expression for the curls is nontrivial (see Lemma 1 below). Indeed, it is a very powerful 
result, since at first it is not evident that there should be no partial derivatives of Pi in (4.9). Fortunately 
that is the case. In fact, it is a requirement, since the Pi are elements of and their derivatives do not 
exist in general. The formula follows from the following lemma coupled with the fact that [Pi](so, si) is a 
homogeneous polynomial of total order i in sq and si. 

Lemma 1, Let ipi{so, si) G 'P*(so, si) be a homogeneous polynomial of total order i in sq and si , where 
So and si are arbitrary functions of some spatial variable in with N = 2,3. Then 

V X (ffso, si)(soVsi - siVso)^ = (i + 2)V^i(so, si)Vso x Vsi . 
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Proof. Notice that 


V X (sqVsi — siVso) = Vso X Vsi + sqV x Vsi — Vsi x Vsq — siV x Vsq = 2Vso x Vsi , 
due to V X Vso = V x Vsi = 0. Now, consider a monomial SgSi, so that 

V X ^Sgs5(soVsi — siVso)^ = Sq^iV x (sqVsi — siVsq) + V(sgs5) x (sqVsi — siVsq) 

= 2sos5^so X Vsi + + 6sos5~^^'Si) x (sqVsi — siVsq) 

= 2sos5^so X Vsi + asQ^^s^^oVso x Vsi — ftsQsJ^^siVsi x Vsq 
= (2 “h ft “h 6 )sqS^Vso X Vsi, 

where it is used that Vsq x Vsq = Vsi x Vsi = 0. Then observe that any homogeneous polynomial 
fii-so, si) is composed of monomials of the form SqS^ of fixed total order ft + 6 = i. The result follows. □ 

Next, record the following important remark. 

Remark. Let fxo = 1 — fj-i, where /xi is an arbitrary function of some spatial variable in with N = 2,3, 

and where p is the order in the coordinates {po, pi). Then for aZZ f = 0,..., p — 1, 

Ef{po,pi) = Pi{pi)Vpi , V X Ef{po, pi) = 0. (4.10) 

With this new ancillary function in our toolset, the i7(curl) shape functions for edge 12 are written 
analogously to those in and the same logic of project —)■ evaluate —)■ blend applies, 

E!{0 = Po{C2)Ef{poi{Ci)) ■ 



evaluate 

In general, the edge shape functions and their curls are 

E!i() = PcmEfipoiiCa)), V X EfiO = Vp.iCb) X EfipoiiCa)) , (4.11) 

where i = 0,... ,pa — 1, {a,b) = (1, 2), (2,1) and c = 0,1. There are pa functions for each edge, giving a 
total of 2p + 2q edge shape functions. 

4.2.2 H (curl) Face Bubbles 

In general, the idea is to consider the tensor product of the ancillary functions Ef and (/>?, evaluated at 
the ID affine coordinate pairs poi{Ci) and /2oi(^2)- To cover all possibilities, one must consider both 
^]‘{Toi{^ 2 ))Ef{poi{Ci)) and (/Uoi(?i))TZf (/Ioi(6))- Clearly, both of these cases are actually generated 
by a single key operator E^, which is defined generally next. 
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Definition. Let (sq, si) and {to, ti) be two pairs of coordinates which are arbitrary functions of some spatial 
variable in with N = 2,3. Let ps be the order in the (sq; si) coordinates, and pt be the order in the 
{to,ti) coordinates. Then 

E°j{so,si,to,ti) = (j)f{to,h)Ef{so,si), (4.12) 

for i = 0,... ,ps — 1 and j = 2,... ,pt. The curls, understood in are 

V X = fj{to,h)V X £;f(so,si) + V4>f{to,ti) x .Ef(so,si) • (4.13) 

Using (3.7), and proceeding as with the 77(curl) edge functions, it is clear that this ancillary operator, 
when evaluated at jdoi{^ 2 )) or (/roi(^ 2 ),satisfies the necessary vanishing trace at all 

edges. There are two families, which together comprise p{q — 1) + q{p — 1) bubble functions. 

Family I: The shape functions for the first family and their corresponding curls are 

EljiO = Wife)), V X 4.(0 = V X fe°(wife), Wife)), (4.14) 

for f = 0,..., p — 1 and j = 2,... ,q. There are p{q — 1) shape functions in this family. 

Family II: The shape functions for the second family and their corresponding curls are 

4.(0 = i?°(wife), Wife)), V X 4(0 = V X 4(wife), Wife)), (4.15) 

for f = 0,... ,q — 1 and j = 2,... ,p. Notice the only difference with the first family is that the entries 
Wi(Ci) Wi(^ 2 ) are permuted (along with the associated orders p and q). Hence, there are q{p — 1) 
shape functions in this family. 


4.3 i?(div) Shape Functions 

By the definition of the H{div) space (see (1.4)) in two dimensions, it is clear that it is isomorphic to 
77(curl). Indeed, the i7(div) shape functions will be the rotation of the corresponding i7(curl) shape 
functions. More explicitly, given a shape function E G i7(curl) and its curl V x 77 G L^, the corresponding 
H (div) shape function with its divergence is 



V - V = V X E. 


(4.16) 


Note in this case the original polynomial space ^ for H (curl) simply becomes the 

i7(div) conforming space x 
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4.4 Shape Functions 


As expected, they are the tensor products of the ID shape functions, and there are pq such functions 
spanning 

4.4.1 Face 

The coordinate free L? shape functions for the quadrilateral faces are 

= [^*](m(ei))[Pi](m(6))(V/ri(ei)xV/ii(6)), (4.17) 

for i = 0,... ,p — 1 and j = 0,..., q — 1. There are pq face functions. The factor V/xi(^i) x V/ri(^ 2 ) 
makes the expression coordinate free (with the ^ coordinates it is 1). 

4.5 Orientations 

For 2D quadrilaterals, only edge orientations need to be considered to ensure compatibility. However, in 
3D elements, quadrilateralwill have the notion of orientation, and one will need to consider this to 
ensure full compatibility. In view of this, first we will show how the concept of edge orientations, already 
explained in §3.3.1, applies to the quadrilateral element. We assume that section has been covered. After 
that, the quadrilateral face orientations are introduced as a preview to the 3D elements. 

4.5.1 Edge Orientations 


‘6 

[0] - - - -> 

3 


a 

A 


m 

[0] - - - -> m 

2 




Figure 4.3; Master quadrilateral with local edge orientations. 

The master quadrilateral has a predefined local orientation for each edge, which reperesents the o = 0 
case. These are illustrated in Figure 4.3. They are our choices for the local orientations. 

Each edge has a local orientation described by the fixed ordering [[]. This induces a master element 
local edge vertex-ordering, which in turn determines a locally ordered pair of affine coordinates, since, over 
a given edge, each master element vertex is linked to one affine coordinate. For instance, over edge 12, 
vi is linked to /ro(?i) and V 2 linked to /ri(^i). On this edge the induced master element local ordering 
is vi-^V 2 (see Figure 4.3), meaning that the locally ordered pair is /7oi(?i) = (fo('^i))F i(Ci))- The 
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locally ordered pair represents the local coordinates, and it serves as the input of the edge local-to-global 
transformation which transforms them to a globally ordered pair (depending on the parameter o). The 
globally ordered pair is then introduced into the edge ancillary functions of the edge shape functions, and 
the resulting functions are said to be orientation embedded shape functions. Thus, for example for edge 12 
the orientation embedded edge functions are. 


(0 = do{i2)<t>f o o-^(/2oi(6)) = 


(o-Q(m(?i))) = M(6)</'f(/^o(6),/^i(6)) ifo = 0; 

= (Aii('^i),fio( 6 )) ifo = l. 




for i = 2,... ,p. This composition with naturally applies to all 2D edge functions in and Tf(curl) 
and their differential forms. Hence, should be composed with cpf, V(/)f, Ef and V x Ef in (4.4) and 
(4.11). 



78- 


-29 




In local edge coordinate 

Figure 4.4: A mesh with global edge orientations followed by transformations to master element and local edge 
coordinates. 


Now, to see a more explicit example starting from the global mesh, consider Figure 4.4. There, one 
can observe a quadrilateral in a mesh with vertices 11, 29, 51 and 78. In 2D, Szabo’s approach (implicitly) 
involves specifying a face vertex-ordering at the mesh that determines the mapping to the master element, 
which has a fixed master element ordering m —?;2 —f 3 —> r’ 4 . However, there is no independent edge vertex¬ 
ordering of each edge. In this case, the face vertex-ordering at the mesh is shown to be 78 —)• 29 —)• 51 —)• 11, 
and at the master element level it coincides with ui —^2 —f 3 —> ^ 4 . The novelty here is that additionally the 
mesh edge with vertices 78 and 29 also has a global orientation given by the edge vertex-ordering 29—)■ 78. 
This edge is mapped to the master element edge 12 (as induced by the face vertex-ordering), and receives 
the induced master element global ordering V 2 ^ vi. Clearly, the local orientation, given by the induced 
local ordering Vi-^V 2 , does not coincide with the global orientation, meaning that the orientation parameter 
is o = 1 for this master element edge. Therefore, the orientation embedded edge shape functions would 
specifically be, 

=Fo(6)</>f ocro(F0i(6)) = Fo(6)<?^f(Fi(6),Fo(6)), 
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for i = 2,... ,p. The neighboring element in the mesh (also sharing vertiees 29 and 78) might have a 
different faee vertex-ordering, but it has the same edge vertex-ordering at the mesh (the global orientation). 
This ean result in another orientation parameter at the master element edge of the neighboring element, but 
by eonstruetion, the faet remains that when mapped baek to the global mesh, the shape funetions will be 
fully eompatible along that shared edge. 

4.5.2 Quadrilateral Face Orientations Explained 



At the mesh In global face coordinates 


In local face coordinates 



Figure 4.5: Quadrilateral face orientations. 

In 3D, as with edge orientations, eaeh faee at the mesh must be given its own global face orientation 
to ensure full eompatibility aeross the boundaries. For quadrilaterals, this is represented by the global 
quadrilateral faee eoordinates or equivalently, by the global face vertex-ordering. For 
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example, given a quadrilateral faee in the mesh, a vertex-ordering of the form a ^ b ^ c ^ d means 
the origin of is located at a, points from a to b, and H 2 points from a to d. Meanwhile, at the 
master element, the mapped face has its own fixed local orientation. It is represented by the coordinates 
or equivalently by the fixed local ordering of fhe form In general, fhe 

fwo sysfems of coordinafes will nol mafch, and Ibis mismafch is represenfed by fhe orienfafion parameter 
o. In facl, fhere are eight possible orienfafions for quadrilaferal faces, meaning o = 0,..., 7. These are all 
illusfraled in Figure 4.5. 

As wifh edges, fhe idea is fo have a local-fo-global Iransformalion which depends on fhe orienfafion 
parameter o. Here, fhe local orienfafion is represenfed by fhe locally ordered quadruple {sq, si,to,ti) 
composed of fhe fwo pairs (sq, si) and (fo, fi)- Thepair, (sq, si), is a locally ordered pair corresponding 
fo fhe^zr^t local face coordinafe (which as an edge coordinafe has fhe local ordering or5]--)-[3]). 

The second pair, (fo,fi), is a locally ordered pair corresponding fo fhe second local face coordinafe 
(which as an edge coordinafe has fhe local ordering |T|-^ [3 or Meanwhile, fhe global orienfafion 

is analogously represenfed by a globally ordered quadruple composed of fwo pairs. The firsf pair is a 
globally ordered pair corresponding fo fhe firsl global face coordinate Sp, and similarly wifh fhe second 
pair. For example, looking af Figure 4.5, when o = 1, SP is associated fo (to, ft), while SP is associated fo 
(si, So), so fhaffhe globally ordered quadruple is (fo, ft, si, sq ). This way, fhe local-fo-global Iransformalion 
is aclually a permulalion dependenl on o lhal can easily be defermined by looking af Figure 4.5. 


Definition. Let so, si, fo and fi be arbitrary variables, and let o = 0,1, 2, 3,4, 5, 6, 7 be the quadrilateral 
face orientation parameter. The quadrilateral face orientation permutation function, a^, is defined as 


'^^('SO, Sl, fo, ft) 


O'?(^ 0 , Sl, fo, ft) = (so, Si, fo, ft) 

0 

II 

0 

C’'p(so, Si, fo, fl) = (fo, fl. Si, So) 

if 0 = 1 

0 'P(so,Sl,fo,fl) = (si,So,fl,fo) 

if 0 = 2 

Clp(so. Si, fo, fl) = (fl, fo. So, Si) 

0 

II 

CO 

O'P(so. Si, fo, fl) = (fo, fl. So, Si) 

if 0 = A 

0'P(so,Sl,fo,fl) = (si,So,fo,fl) 

if 0 = 5 

Clp(^ 0 , Si, fo, fl) = (fl, fo. Si, So) 

II 

0 

^O-p(so. Si, fo, fl) = (so. Si, fl, fo) 

0 

II 


(4.18) 


As wifh edges, all lhal is required is fo compose fhe quadrilateral face ancillary functions and Iheir 
differenlial form (Ihose wifh superscripl □) wifh fhe local-fo-global Iransformalion given by a^. This should 
be done in all 3D shape funclions associaled fo quadrilateral faces. More concrele examples will be given 
in fhe 3D elemenfs as fhe documenl progresses. 
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5 Triangle 



Figure 5.1: Master triangle with numbered vertices. 

The triangle is the 2D simplex. The master element for triangles in the x = (xi, X 2 ) space is the set 
{x G : xi > 0, X 2 > 0, xi + X 2 < 1}- It is illustrated in Figure 5.1. 

Denote vertex a by Va, so that vq = (0, 0), vi = (1,0) and V 2 = (0,1). As described in §1.9, the 2D 
affine coordinates, vq, ui and 1 ^ 2 , can be easily calculated for this master triangle: 

^^{x) = 1 — Xl — X2 , l'l{x) = Xl , V2{x)=X2. (5.1) 

Their gradients are 

Vt^o(x) = Vz.i(x) = (j), Vv 2 {x) = (l). (5.2) 

Like the quadrilateral and segment, the triangle exhibits a correspondence of its vertices and its affine 
coordinafes. Looking af fhe formula x = vq{x)vo + z^i(x)xi + i' 2 {x)v 2 this relation is evident, in the sense 
that each vertex is linked to one affine coordinafe (ifs corresponding weighf). For example xi is linked fo 
the affine coordinafe i^i(x), and indeed if fakes fhe value 1 when x = vi while if is zero af fhe ofher fwo 
vertices. 

Exact Sequence 

As wifh fhe quadrilaferal, fhe friangle will have 2D discrefe polynomial exacf sequences fhaf represenf fhe 
confinuous exacf sequence (1.2) and ifs rofafed analogue (1.5). They are 

pp j\[p 

’ (5.3) 

-pp ^ 

where = V^{xi,X 2 ) is fhe space of polynomials of fofal order p. The spaces and TZT^ are fhe 
Nedelec and Raviarf-Thomas spaces for simplices: 

j^v = (pp-i)Af © |p e {j5p)N . ^ = 0 for all x G M^} , (5.4) 

pjpp ^ (pP-i)Af © |p e (pP)N . with e pp-1 and X G M^} . (5.5) 
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In the case of the triangle, the number of spatial dimensions is = 2. Note the sequence has an overall drop 
in polynomial order of one. This makes it compatible with the construction presented for the quadrilateral. 
Moreover, all of the spaces in the exact sequences above are invariant under affine transformations. This 
implies the exact sequence takes the same form for any given triangle (provided it is mapped via an affine 
fransformafion from fhe masfer friangle, which is always possible). 

5.1 Shape Functions 

All shape funclions defined here will lie in V^, which has dimension + 2)(p + 1). Moreover, a careful 
counf of fhe linearly independenf shape funclions will coincide wilh lhal dimension, so lhal indeed fhe space 
is spanned. 

5.1.1 Vertices 

In Ibis case, for a given vertex, fhe associaled shape funclion will be ils related affine coordinale. For 
example, vq is linked lo vo{x), so ils associaled verlex funclion is simply 

= iyo{x) . 

As expecled, if vanishes al fhe disjoinl opposile edge 12 and ils Irace over Ihe adjacenl edges is a ID 
vertex funclion associated lo Ihe vertex. Indeed, Ihe Iraces are explicilly, 

yo{x)\x2=Q = 1- Xl = ^oixi) , Vo{x)\l-xi-x2=0 = 0 , vo{x)\xi=o = 1- X2 = noix2) ■ 

Laslly, Ihe funclion decays linearly and is in Ihe lowesl order possible space, so lhal if respecls Ihe 
hierarchy in p. 

In general, Ihe vertex funclions and Iheir gradienl are, 

= I2a{x) , V(t)"'{x) = VVa{x) , ( 5 . 6 ) 

for a = 0,1, 2. There are a lolal of 3 vertex funclions (one for each vertex). 

5.1.2 Edges 

For Ihe conslruclion of Ihe edge funclions consider edge 01 as an example. To salisfy compalibilily wilh Ihe 
quadrilateral, Ihe shape funclions should have {flQi{xi)) as Ihe Irace over edge 01, and be zero al Ihe Iwo 
olher edges. Addilionally, Ihey should be polynomial (Ihey should be in V^). 

Unforlunalely, al firsl sighl Ihe conslruclion is nol Irivial since many inluilive approaches lead lo violal- 
ing some of Ihe required properlies mentioned above. This is in large pari due lo Ihe facl lhal Ihe Iriangle 
is not a Cartesian producl of lower dimensional elemenls. Luckily, Ihese issues are all solved by Ihe pro¬ 
cess of homogenization, which is viewed as a particular extension of Ihe lower dimensional edge functions 
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(ff {11 qi{x\)). The idea is to exploit the 2D triangle affine coordinates associated to edge 01, which are vq 
and ui (those linked to the vertices vq and vi which compose the edge). Indeed, let the shape functions for 
this edge be 


4 >i(x) = (/>f(m(x)) = [Li]{pQi{x)) 


(uoix) + i'i{x)y[Li] 


l^oix) 1^1 (x) \ 

uo(x)+i/i(x) ’ I^o(x)+I^l(x) J ’ 


with i = 2,... ,p. Here, property (2.32) was used. Notice that 1/2 = 0 is the equation for edge 01 and 
similarly with the other edges. Hence, using the vanishing properties of cpf in (3.7), it follows that the trace 
properties are satisfied: 


<^-(x)U2=o = 4>f{i^oi{x))\u2=o = yfimixi)) , 

yi{x)\l-xi-X2=0 = (/>f ( 1 ^ 01 ( 2 :)) 11.0=0 = </>f( 0 ,/il(xi)) = 0 , 
<^|(x)Ui=0 = (t)f{vQi{x))\y^=o = (t)f{po{x2),Q) = 0. 


In this case, the decay of each shape function is represented by the nonlinear function {^^{x) + ui{x)y, 
which comes hidden within the homogenization. Additionally, (^f(i7oi) £ by homogenization 

while vq { x ), vi{x) G V^{x), meaning that the edge functions (j)y{vQi[x)) are in 'P*(x) C 'PP(x) as required. 




Figure 5.2: Edge projection from P to P', and the logic project —evaluate —> blend. 

As with the quadrilateral, there is a geometrical interpretation to these expressions, and again it follows 
the fundamental logic of project —)■ evaluate —blend, which is marked below: 


yt{x) = (^7oi(x)) = fro(x) + Z^l(x))>f ( , U,{x)+Mx) 


I'ljx) 


blend 


project 


evaluate 


The coordinates (projected ID coordinates, because they sum to 1 for all x. Note this is not 
true for the coordinates (z/q, z^i) in general (only over the edge itself). As argued before for the quadrilateral, 
the projected coordinates represent a point over edge 01: 


(xi,X2) ^ (t^>0). 

Geometrically, it consists of finding the intersection P' = (pz^,0) of the edge with the projecting line 
passing through the original point P = {xi,X 2 ) and the disjoint opposite vertex. This projection and the 
logic of the construction is illustrated in Figure 5.2. 
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The complete list of edge functions with their gradients is 

(Piix) = (pfil^abix)) , = V(/)f (z7afe(x)) , (5.7) 

with i = 2,... ,p, and 0 < a < 6 < 2 (so (a, 6) = (0,1), (1, 2), (0, 2)). As usual, there are a total of p — 1 
edge functions for every edge, for a total of 3(p — 1) edge functions. 

5,1.3 Face Bubbles 

These intuitively involve all three affine coordinates {uq, ui, z^ 2 )-The first natural idea is to use the existing 
edge shape functions and multiply by another function which vanishes at the remaining edge. That is, 

(f^ijix) = (l)f{iyo{x),iyi{x))^j{u2{x)) . 

Here, (j)j{y) should also be an function with domain y G [0,1] (since 0 < V 2 {x) < 1) and which 
vanishes at y = 0, so that i^j(O) = 0. Notice that in the quadrilateral, when constructing i;/)?- the choice 
was naturally (pj = Lj since it was required that <?!)j(0) = = 0, so it should vanish at both endpoints. 

For the triangle there is much more liberty for the choice of (pj, and indeed it can be chosen such that it has 
many advantages. Following Beuchler and Schdberl (2006) it is chosen as a Jacobi polynomial, (pj = L“, 
with a = 2i, where i is the order of (p^. This results in 

4>\j{x) = (pf{oo{x),vi{x))Lf{u 2 {x)). 

Thinking ahead to the tetrahedron, this definition can be easily generalized. Just as with (pf, the gener¬ 
alization is nothing more than the homogenization of the previous formula, as written in (2.37). 

Definition. Let sq. si and S 2 be arbitrary functions of some spatial variable in with N = 2,3, and 
denote by ps the order in the coordinate triplet (sq, si, •52)■ Then 

4>ijiso, Si, S 2 ) = <pf{so, si)[Lf]{so + Si, S 2 ) = [Li,Lf]{so, si, S 2 ), 

for i > 2, j > 1 and n = i + j = 3,... ,ps. The gradients, understood in are 

V(pfj{so, Sl,S 2 ) = [Lf]{so + Si, S 2 )V(^f (so, Si) + (pf{so, si)V[L|*](so + Si, S 2 ), 

where for any a, 

K'^O + 'Sl, S 2 ) = [Pj^]^](so + Si, S 2 )VS 2 + [i?“_]^](so + Si, S 2 )V(S 0 + S 1 + S 2 ) 

= P^-i{s2] So + Si+S 2 )Vs 2 + R^_i{s 2 ', Sq+ Si + S 2 )V(so + Si+S 2 ) . 


(5.8) 

(5.9) 


(5.10) 


Note the indexing was shown as i > 2, j > 1 and n = f + y = 3,... ,ps. In many hp codes it is useful 
to enforce hierarchy, so that the shape functions are organized by total order. That is, first list all order 3 
bubbles, then all order 4 bubbles, and so on. Hence when coding is in mind, it is useful to have an outer 
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loop with numbering n = 3,... ,ps and an inner loop indexing with either i = 2,... ,n — 1 (so j = n — i) 
or j = 1,..., n — 2 (so i = n — j). 

Regarding the vanishing properties of this new aneillary funetion, it suffiees to rewrite (2.40) so that for 
any sq; si, S 2 , and all i > 2, j > 1, 

<i>tjiso,si,0) = (/)g(so,0,S2) = (j)tjiO,si,S 2 ) = 0. (5.11) 

As expeeted, when the eoordinates are 2D affine eoordinates, the formulas for (j)fj and its gradient are 
simplified. For fhis, reeord fhe nexf remark. 

Remark. Let vq = 1 — vi — 1 / 2 , where vi and 1^2 tire arbitrary functions of some spatial variable in 
N = 2, 3, and where p is the order in the coordinates (uq, vi, 02 ). Then for all i > 2, j > 1, and 
n = i + j = 3,...,p, 

01 , 02 ) = 4>f{oo, 0 i)Lf{02) , 

01 , 02 ) = Lf(o2)Vff{oo,oi) + 4>f(oo,oi)P^^{o2)Vo2. 

In general, fhe friangle faee bubbles and fheir gradienfs are 

(l^ljix) = </>§(F'oi 2 (x)) , V(/)^j(x) = V(/>g(Foi 2 (x)), (5.13) 

where i > 2, j > 1 and n = i + j = 3,... ,p. The vanishing eondifions are salisfied by eonsfruefion or 
simply by looking af (5.11). There are a fofal of — l)(p — 2) faee bubbles for fhe friangle. 

5.2 i?(curl) Shape Functions 

The dimension ofM^ in fwo dimensions is p(p + 2). A eareful eounf of fhe linearly independenf eonforming 
shape funelions fo be presenfed fhroughouf Ibis seefion will eoineide wifh fhaf dimension. Showing fhaf fhe 
funelions eonsfruefed are in is nonfrivial, buf will follow from fhe nexf insfrumenfal lemma. 

Lemma 2. Let x G for N = 2,3, and fn G V^(x) be any polynomial of total order n in the coordinates 
X = (xi,..., xtv)- Given sq and si coordinates in 'P^(x) (linear functions in x), it follows that the Nedelec 
space of order n + 1, contains the function 

fn{-){soVsi-SiVso) GAA"+1. 

Proof. Reeall fhe definilion of fhe Nedelee spaee, 

J^P ^ (pp-l)Af 0 ^ ^ Q ^ ^ ^iVj _ 

The (affine) eoordinafes are linear funelions in x = (xi,..., xtv), so fhaf 

Sk(x) = ak + bk-x, 


41 



for Ofc G M, G and /c = 0,1. As a result Vsfc(x) = bk, and 

so(x)Vsi(x) - si(x)Vso(a;) = (ao^i - ai&o) + ((&o • x)bi - (6i • x)bo) . 

" -V-^ '-V-^ 

=A =B(x) 

Clearly, A G while B G (P^)^. Moreover, 

X ■ B(x) = X ■ ((bo ■ x)bi — (6i • x)bo) = (bg ■ x)(bi ■ x) — (bi ■ x)(bo ■ x) = 0, 

so that B G {E G (P^)^ : x ■ E(x) = 0}, and sqVsi — siVso £ 

Now, fn G P'^ = P^~^ 0 "P”, for n > 1 ean always be deeoupled into /„ = fn-i + fn, where 
fn-i G and /„ G P”. As a result 

fn(x)(so(x)\/Si(x) - Si(x)Vso(x)) = fn(x)A + fn-l(x)B(x) + fn(x)B(x) , 
where it is clear f^A + fn-iB G (P^)^ and /„P G (p^+i)-^. Meanwhile, 

X ■ (fn(x)B(x)) = Jn(x)x ■ B(x) = 0 . 

Hence,/„(• )(soVsi — siVso) £ □ 


5.2.1 H (curl) Edges 


Having seen what ocurred with edge functions due to the process of homogenization, it is wise to 
consider the general definition of P(curl) edge functions in (4.8), using 2D affine coordinafes as entries. 
Indeed this will suffice. For instance, for edge 01, the shape functions are 


P®(x) = Ef(voi(x)) = [Pi](voi(x)){vo(x)yyi(x) - vi(x)Vvo(x)^ 

)pf(Foi(x)) 


= (i'o(x) + ui(x)y[Pi] 


l'o(x) 


i'i(x) 


i/o(x)+ui(x) ’ uo(x)+iyi(x) 


= (l^o(x) + l^i(x)r+]Ef{ uo(x)+Z(x) 


l'o(x) 


‘'l(x) 


blend 


project 


evaluate 


for i = 0,... ,p — 1. By Lemma 2 it easily follows P| G C AfP as desired. The tangential trace 

properties are also satisfied. To see this, it suffices to look at Eq, since Ef(x) = [Pj](Foi(x))PQ(x). Its 
traces are 


tr(P^(x))U 2 =o = E^(uoi(x))y 2 =o • (^1 - = Vui(x) ■ (xi - xq) = 1 , 

tr(Po(a;))|l-xi-a;2=0 = EQ(i 7 oi(x))\uo =0 ■ (V2 - Ui) = -ui(x)Viyo(x) ■ (V2 - Xi) = 0 , 

tr(Po(a:))Ui=0 = EQ(i2oi(x))\ui=0 ■ (V2 - Vg) = l^o(x)Vui(x) ■ (V2 - Vg) = 0 . 
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Here the vanishing traees follow from the faet that the gradient of a funetion is always perpendieular to the 
tangents of its isosurfaees. Henee is perpendieular to the tangents of the set where i'o{x) = 0, whieh 

is preeisely where V 2 — vi lies. Regarding the nonzero traee, the formula Vui{x) ■ (ui — vq) follows after 
replaeing uq = 1 — ui, and noting that V{iyo{x) + viix)) ■ (ui — uq) = —Vv 2 {x) ■ (ui — vq) = 0. 

More generally, the edge funetions and their eurls are 

Ff(x) = (i7a6(x)), V X B^(x) = V X Ef(,7,Ux)), (5.14) 

with i = 0,... ,p — 1, and 0 < a < 6 < 2 (so (a, 6) = (0,1), (1, 2), (0, 2)). There are a total of p edge 
funetions for every given edge, giving a total of 3p edge funetions. 

5.2.2 H (curl) Face Bubbles 

The eonstruetion of Tf(curl) faee bubbles is parallel to that of bubbles. The idea is to multiply Tf(curl) 
edge funetions with an funetion vanishing at zero. Again, it is ehosen aeeording to Beuehler et al. (2013) 
as the Jaeobi polynomial, L“, with a = 2i + l, where i is the order of Ef. Like the quadrilateral, the triangle 
has two elosely related families generated by the same aneillary operator defined below. 

Definition. Let sq. ■si cind S 2 be arbitrary functions of some spatial variable in with N = 2,3, and 
denote by ps the order in the coordinate triplet (sq, si, •52)■ Then 

E^^{sq, Si, S 2 ) = [Lf^\so + si, S 2 )£'f (so, si), (5.15) 

for ^ A 0, j > 1 and n = i + j = 1,... ,ps — 1. The curls, understood in are 

Vx.eA(so, si, S 2 ) = [Lf+^](so+si, S 2 )V xF;f (so, si)+V[Lf+^](so + si, S 2 ) xT;f (so, si), (5.16) 

where V[L^*^^](so + si, S 2 ) is computed from (5.10). 

As with the edge funetions, it is elear that when evaluated in any permutation of (i/q, , 1 ^ 2 ) the vanishing 

traee properties are satisfied. The fwo families defined below give a fofal of p{p — 1) Lf(curl) friangle faee 
bubbles (eaeh family has ^p{p — 1)). 

Family I: The shape funetions for fhe firsl family and fheir eurl are 

Ejj{x) = Ef-{Poi2 {x)), V X Ejj{x) = V X Ef-^douix )), (5.17) 

for i > 0, y > 1 and n = i + j = 1,... — 1. By Lemma 2, ET G where n = i + j. There are 

^p{p — 1 ) bubble funetions in Ibis family. 


43 



(5.18) 


Family II: The shape funetions for the seeond family and their eurl are 


eUx) = E^A Pi2o{x)) 


V X eUx) = V X E^Avi2o{x )), 


for i > 0, j > 1 and n = i + j = l,...,p — 1. The only differenee with the first family is that the entries 
are permuted to i7i2o(a;) = {vi{x), V 2 {x), uo{x)) of vqi 2 {x) = {iyo{x),i'i{x),i' 2 {x)). Again, they 

lie in Efj G where n = i + j, and there are ^p{p — 1) bubble funetions in this family. 

5.3 i?(div) Shape Functions 

In two dimensions, the Tf(div) spaee (see (1.4)) is isomorphie to H{curl). Indeed, just like with the quadri¬ 
lateral, given a shape funetion E G H{cur\) and its eurl V x i? G L^, the eorresponding iT(div) shape 
funetion with its divergenee is 



V - v = V X E. 


(5.19) 


Although not immediate, note that in two dimensions the original polynomial spaee for H (curl) simply 
becomes the H (div) conforming space TZT^ after the rotation, as required. 

5.4 Shape Functions 

These span V^~^, so there should be \p{p + 1) linearly independent shape functions. 


5.4.1 

Again, carefully chosen Jacobi polynomials are present in this construction. The shape functions are 

i’ljix) = [Pi]{vo{x),Vi{x))[Pf"^\vQ{x) + Vl{x),V2{x)) 

= [Pi,Pf+\vm2{x)){Vvi{x)xV02{x)) , 

where i > 0, j > 0 and n = i + j = 0,...,p—1. Clearly the functions lie in ^ and there are \p{p + f) 
such functions. The factor Vv'i{x) x Vr' 2 (x) makes the expression coordinate free (with the x coordinates 
it is 1). 


(5.20) 


5.5 Orientations 

Only edge orientations are required to ensure compatibility of 2D triangles. However, in 3D elements, 
triangle/ace^ have the notion of orientation, and one needs to consider this to ensure full compatibility. 
In §4.5.1 it was shown how to apply the principles introduced in §3.3.1 to construct orientation embedded 
edge functions. It is convenient to have read those sections. In what follows, a brief illustration of how it 
analogously applies to the triangle is given in §5.5.1. Afterwards, in §5.5.2, the triangle face orientations are 
described. 
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5.5.1 Edge Orientations 



Figure 5.3: Master triangle with local edge orientations. 

The master triangle has a predefined local orientation for each edge, which reperesents the o = 0 case. 
These are shown in Figure 5.3. They are our choices for the local orientations. 

As with the quadrilateral, the local orientations induce a master element local edge vertex-ordering 
which then determines a locally ordered pair of affine coordinates. In the triangle, the correspondence 
between vertices and affine coordinates is trivial, since all one needs to know is that the vertex Va is linked 
to the affine coordinate Ua, where a = 0,1, 2. As an example take the local edge orientations shown in Figure 
5.3. Then, it is clear that the induced local edge vertex-orderings we vq-^vi,vi-^V 2 and vq-^V 2 for edges 
01, 12 and 02 respectively. Therefore, the locally ordered pair of coordinates are {vq,vi), {ui,V 2 ) and 
(t'O) t^ 2 )- These are then transformed to a globally ordered pair via the edge local-to-global transformation 
do , and inputted into the different ancillary functions to give the orientation embedded shape functions. For 
example, for edge 01, the orientation embedded edge functions are 


oa^{voi{x)) 


</’f (^0 (^^oi(a:))) = (/)f(i7oi(x)) if o = 0 , 
= (/)f (Fio(x)) if o = 1, 


for t = 2,... ,p. For a more detailed example which analogously applies to the triangle see §4.5.1. 


5.5.2 Triangle Face Orientations Explained 

In 3D, each face at the mesh must be given its own global face orientation to ensure full compatibility across 
the boundaries. For triangles, this is represented by the global triangle face coordinates = (H^, H^), or 
equivalently, by the global face vertex-ordering. For instance, given a triangle face in the mesh, a vertex¬ 
ordering of the form a ^ b —?■ c means the origin of is located at a, points from a to b, and 
points from a to c. Meanwhile, at the master element, the mapped face has its own fixed local orientation. 
It is represented by the coordinates or equivalently by the fixed local ordering of the form 

In general, the two systems of coordinates will not match, and this mismatch is represented 
by the orientation parameter o. For triangle faces, there are six possible orientations, meaning o = 0,..., 5. 
These are all illustrated in Figure 5.4. 
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Figure 5.4: Triangle face orientations. 

Like quadrilaterals and edges, one needs a loeal-to-global transformation dependent on o. For this, the 
local orientation is represented hy the locally ordered triplet (sq, si, •52)- The global orientation is analo¬ 
gously represented hy a globally ordered triplet naturally indueed from the global face vertex-ordering. For 
example, looking at Figure 5.4, when o = 1, the global ordering a ^ b ^ c eorresponds to the ordering 
Q] — )■ [2] — )■ [0], whieh in turn eoresponds to the globally ordered triplet (si, S2, sq)- Henee, onee again the 
loeal-to-global transformation ean be represented by a permutation, a^, dependent on o. It ean be deter¬ 
mined by observing Figure 5.4. Subsequently, all that is required is to eompose the triangle face aneillary 
funetions and their differential form (those with superseript A) with the loeal-to-global transformation 
This should be done in all 3D shape funetions assoeiated to triangle faees. More eonerete examples will be 
given later. Finally, the permutation funetion, a^, is preeisely defined next. 
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Definition. Let sq. si S 2 be arbitrary variables, and let o = 0,1, 2, 3,4, 5 be the triangle face orienta¬ 
tion parameter. The triangle face orientation permutation function, a^, is defined as 


O'o ('So, 'Si, S2) 


(T^(S0,'S1,S2) = (S0,S1,S2) 
O-f (so,'Sl,'S2) = (S1,S2,'S0) 
<7^(S0,'Si, S 2 ) = (•S2,-So,'Si) 
£7^(S0,'S1,S2) = ('S0,S2,'Sl) 
O-f (S0,'S1,S2) = ('Sl,S0,'S2) 
'Si, S 2 ) = ('S2,'Sl,So) 


if 0 = 0, 
if 0 = 1, 
if 0 = 2, 
if 0 = 2,, 
if 0 = 4:, 
if 0 = 5. 


(5.21) 
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6 Hexahedron 



Figure 6.1: Master hexahedron with numbered vertices. 

The master element for hexahedra is (0,1)^. It is shown in Figure 6.1 in the ^ = (^i, ^ 2 , ^s) space. The 
master hexahedron is the Cartesian product of three segments. 

There are three pairs of ID affine coordinates: 

fo(6) = i-6) fi( 6) = 6 =► V/io(Ci) = ^ 0 ^ 

/^o(6) = i-6, fii(6) = 6 =► V/io(6) = , V/xi(^2) 

Fo(6 ) = 1-6) ^1(6)= '^3 => V/io(?3) = ^ 0^^ , V/ri(^3) 

These will be used explicitly or implicitly in the formulas that follow. 

Just as with quadrilaterals, there are natural relationships between vertices, edges and faces, and the 
affine coordinates. In fact, each vertex is linked to three affine coordinates, each edge is linked to two affine 
coordinates, and each face is linked to one affine coordinate. The linked affine coordinates take the value 1 
at the associated topological entity. For example, vertex 1, vi = (0,0,0), is linked to the affine coordinates 
Fo(?i)! fto(^ 2 ) and /ro(^ 3 )> edge 12 is linked to the affine coordinates /Uo(C 2 ) and fioiCs), and face 1234 is 
linked to affine coordinate 



Exact Sequence 


Recall the 3D exact sequence for simply connected domains (1.7). The corresponding discrete polynomial 
exact sequence is of the form 




Vx 


QP,q,r '' yp,q 


YP’1’‘ 


(6.2) 
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where the standard Nedelec’s spaces Nedelec (1980) of the first type for the hexahedron are utilized: 

WP,q,r ^ Qp,q,r ^ j,p^ ^ ^ 

QP,Q,r _ Qp-l,q,r ^ Qp,q-l,r ^ Qp,q,r-1 ^ 
yp,q,r _ Qp,q-l,r-l ^ Qp-l,q,r-l ^ Qp-l,q-l,r ^ 


(6.3) 


As with the quadrilateral, there is a natural anisotropy of the element, which has order p in the direction, 
q in the ^2 direction and r in the ^3 direction. The hierarchy should be maintained in p, q and r separately. 
It will sometimes be convenient to refer to pa as the order in the direction, so that pi = p, P 2 = q and 
P3 = r. 

6.1 Shape Functions 

It will be clear that the {p + l)(q + l)(r + 1 ) shape functions defined in fhis secfion lie in Qf’'?’’’ and span 
fhe space. 

The ideas in fhis secfion are fhe same as wifh fhe quadrilaferal (see §4) buf in fhree dimensions. This 
will simply franslafe fo adding an exfra blending funclion fo accounf for fhe exfra dimension. Hence, fhe 
frace properties will nol be analyzed in defail as fhey easily follow. 


6.1.1 Vertices 

Wifhouf any delays, fhe vertex shape funclions and fheir gradienf are 

= pa{(l)Pb{C2)fJ‘c{^3) , 

for a = 0,l, 6 = 0,1 and c = 0,1. There are a fofal of 8 verfex funclions (one for each verlex). 


(6.4) 


6.1.2 Edges 

Again, fhis is analogous fo fhe quadrilateral case, buf wifh an exfra blending function. Take for example 
edge 12. Then, fhe shape funclions are 

= mo(6)mo( 6) y>f , 

project 


blend 


evaluate 


for f = 2 ,..., p. The projection being implied is: 

( 6 , 6 ;? 3 ) I—> (6)6)0) I—> (6,0,0) 
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It is illustrated in Figure 6.2. It consists simply of finding the intersection P" = (.^i, 0,0) of the edge with 
the normal plane passing through the original point P = (^i, ^ 2 , ^s)- Alternatively it can be interpreted in 
two steps. First it is projected to the closest point P' in an adjacent face, by using the normal to the face. 
Once in the face, it is projected again to the desired edge using the traditional two dimensional quadrilateral 
projection (see Figure 4.2). 




Figure 6.2: Face projection from P to P' followed by an edge projection from P' to P". 
The general edge shape functions and their gradient are 


4’t{€j = /^e(6)/rd(6)(/'f (FOl(Ca)) , 

where i = 2,... ,pa, (a, b, c) = (1, 2,3), (2, 3,1), (3,1, 2), d = 0,1 and e = 0,1. For example for edge 12, 
this would correspond to (a, b, c) = (1, 2, 3), d = 0, e = 0 and pa = p, for edge 23 it is (a, b, c) = (2, 3,1), 
d = 0, e = 1 and pa = q, and so on. For each edge there are pa — 1 shape functions, leading to a total of 
4(p — 1) + 4(q — 1) + 4(r — 1) edge functions. 


6.1.3 Faces 

Again, by adding as a blending factor the associated affine coordinafe fhe consfrucfion becomes frivial. For 
example, for face 1234, fhe shape functions are 

(C) = 4>?j ( /^oi(6)^Foi(6p , 

blend project 

^-V-^ 

evaluate 

where i = 2,... ,p and j = 2,... ,q. The projection is already illuslrafed in Figure 6.2, and simply consisfs 
of finding fhe infersecfion P' = (■^i, ■^ 2 ) 0) of fhe normal fo fhe face fhaf passes fhrough fhe original poinf 
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In general the faee shape funetions and their gradient are 






(6.6) 


where i = 2,... ,pa, j = 2,... ,p6, (a, b, c) = (1,2, 3), (2,3,1), (3,1, 2), and d = 0,1. For example for 
faee 1234, this would eorrespond to (a, b, c) = (1, 2, 3), d = 0, pa = p and p^ = q, for faee 1265 it is 
(a, b, c) = (3,1, 2), d = 0, Pa = r and ph = p, and so on. For eaeh faee there are {pa — l){pb — 1) shape 
funetions, for a total of 2{p — l)(g — 1) + 2{q — l)(r — 1) + 2(r — l)(p — 1) faee funetions. 


6.1.4 Interior Bubbles 

These are eonstrueted like the faee funetions, but by using edge bubbles (pf instead of the linear blending 
faetor pd- This will ensure the neeessary vanishing traee properties. 

The interior bubbles and their gradient are 


(6.7) 


= Li{pi{^i))Lj{pi{^2))Lk{pi{^3)) = </>^(/(oi(Ci),/7oi(C2))0f(Foi(6)), 

= </’^(Foi(6),/7oi(6))V(/>f (/2oi(6)) + </>fc(/^oi(6))V(/>°(/2oi(6),/7oi(6)), 
for i = 2,..., p, j = 2,... ,q and k = 2,... ,r. Clearly there will be (p — 1) (g — 1) (r — 1) interior bubbles. 


6.2 i?(curl) Shape Functions 

It will be elear that the p{q + l)(r + 1) + q{r + l)(p + 1) + r{p + l)(g + 1) linearly independent shape 
funetions span x x as required. 

The ideas in this seetion are the same as with the quadrilateral but in three dimensions. This will simply 
translate to adding an extra blending funetion to aeeount for the extra dimension. The strueture of projeeting, 
evaluating and blending still holds in H (curl), and the projeetions (and even the blending funetions) are the 
same as those in As with the analysis of the traee properties will be superfluous. 


6.2.1 H (curl) Edges 

These will just be the quadrilateral edge funetions with the extra blending faetor. They are 

EtiO = d'e{ic)d-d{ib)Ef{floi{^a)) , 

V X Et{i) = [pe{^c)^lid{^b)+^^d{^bW^^e{^c)) X (poi(^a)) , 


(6.8) 


where i = 2,... ,pa, {a, b, c) = (1, 2, 3), (2, 3,1), (3,1, 2), d = 0,1 and e = 0,1. Notiee the form is very 
similar to that of edge funetions. For eaeh edge there are pa shape funetions, giving a total of 4p+4q+4r 
edge funetions. 
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6.2.2 H (curl) Faces 


The pattern goes on, but this time with the two families. There are a grand total of 2{p{q — 1) + q{p — 1)) + 
2 {q{r — 1) + r{q — 1)) + 2 {r{p — 1) + p{r — 1)) faee shape funetions. 


Family I: The shape functions and their curl are 

^iiiO = P'OliCb)) , 


^ X EfjiC) = M<i(^c)V X Polish)) + V^rf(^c) X (/foi(Ca), MOl(6)) 


(6.9) 


where i = 0,... ,pa — 1, j = 2,... ,ph, {a,b,c) = (1, 2, 3), (2, 3,1), (3,1, 2), and d = 0,1. For each face 
there are Pa{pb — 1) shape functions in this family. 


( 6 . 10 ) 


Family II: The shape functions and their curl are 

EijiC) = ld'd{ic)Efj{poi{ib),Uoi{ia)) , 

V X Elj{£) = Pd{ic)y X £'j°(/i0l(^b),F0l(^a)) + '^fJ^diCc) X E'^jipoii^b), PoliCa)), 

where i = 0,... ,pb — 1, j = 2,... ,pa, {a,b,c) = (1, 2, 3), (2, 3,1), (3,1, 2), and d = 0,1. Again, recall 
the only difference with the first family is that the entries {floiiCa) 7 fdoiiCb)) and their corresponding order 
(pa,Pb), are permuted. For each face there are pb{pa — 1) shape functions in this family. 


6.2.3 H (curl) Interior Bubbles 

These can be constructed by using edge bubbles cj)^ as blending functions instead of the linear blending 
Pd in the expressions for the face shape functions. All permutations of (a, 6, c) leading to linearly indepen¬ 
dent functions must be considered. In the end, three famillies corresponding to the cyclic permutations of 
(1, 2,3) comprise the interior bubbles. 

The interior functions and their curl are 

Eijki^) = </>f (Foi(^c))-Ej°(/2oi(^a),m(4)), 

V X = <^f(m(ec))VxF;°(M0i(ea),m(6)) + V</.f(^0l(ec))xS°(M0l(^a),F0l(6)) • 

fori = 0,... ,pa - I, j = 2,...,pb, k = 2,...,pc, and (a, 6, c) = (1, 2, 3), (2, 3,1), (3,1, 2). There will 
be a grand total of p(q — l)(r — 1) + q{r — l)(p — 1) + r{p — 1)((7 — 1) interior bubble functions. 


( 6 . 11 ) 


6.3 i?(div) Shape Functions 

It will be clear that all shape functions lie in x x which has 

dimension rq{p + 1) + pr{q + 1) + qp{r + 1). 

This is the first time that the space dd(div) is tackled in 3D. As expected, it requires of some analysis to 
develop the correct structure at first, but afterwards one can proceed very similarly as the previous spaces. 
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6.3.1 H{dW) Faces 

First, recall from § 1.6 that the normal trace of the H (div) face functions should be a 2D face function. For 
the purposes of motivation, take for instance face 1234. From (4.17), the 2D face functions are of the form 
[-Pi](Foi(^i))[^j](Foi(C 2 ))- Meanwhile, the normal vector to face 1234 is (0, 0,1) = V;Ui(.^i) x V;Ui(.^ 2 )- 
When coupled with a blending factor, ^ 0 (^ 3 )? representing a linear decay (like that of this suggests, 

= Fofe)[^’d(Moi(6))[^,](m(6)) Vmi( 6) X V/ri(6) = Fote)i5^f(m(ei)) X Efimia)), 

for z = 0,... ,p — 1 and j = 0,..., g — 1. In fact, this expression makes a lot of sense, since a function 
normal to the face should be perpendicular to the two tangential H (curl) edge functions. The cross product 
then seems like a natural idea. Indeed, one can laboriously check that the desired normal trace properties 
are satified at all faces. This motivates the definition of a new ancillary operator presented next. 

Definition. Let (sq, si) and {to, ti) be two pairs of coordinates which are arbitrary functions of some spatial 
variable in with N = 3. Letps be the order in the (sq, si) coordinates, andpt be the order in the {to, fi) 
coordinates. Then 

V^{so,si,to,ti) = Ef{so,si) X Ef{to,ti), (6.12) 

for i = 0,... ,ps — 1 and j = 0,... ,pt — 1. Their divergence, understood in is 

V • ■Fj°(so,si,fo,fi) = Ef{to,ti) ■ (V X Ef{so,si)) - Ef{so,si) ■ (V x Ef{to,ti)). (6.13) 


Record also the next useful remark. 


Remark. Let = 1 — and = 1 — p^i \ where and p^^^ are arbitrary functions of some 
spatial variable in with N = 3, and where p(o) are the orders in the coordinates {pg^\ p[^^) 

and {p'^\p^i^) respectively. Then for all i = 0,... ,p(o) “ and j = 0,... ,p(i) — 1, 


^■V^{pf, 


tT: 

tT: 


ho'’ 

ho^ 


rf^) = 0. 


(6.14) 


Finally, the shape functions and their divergence are 

= hd{^c)y^{P0l{ia),P0l{ib)) , 

V • = ^hd{Q • v^{poi{ia),m{ib )), 

where f = 0,... ,pa — 1. J = 0,... — 1, (a, 6, c) = (1,2, 3), (2,3,1), (3,1, 2), and d = 0,1. There are 

PaPb face functions for each face, leading to a total of 2pq + 2qr + 2rp face functions. 
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6.3.2 H (div) Interior Bubbles 


Using the same reasoning as with H (curl) interior bubbles, there will essentially be three families of bubbles 
corresponding to the cyclic permutations of (1,2,3). The interior bubbles and their divergence are 

= </’fc(m(Cc))^i°(m(?a),m(^6)), 

V • u^(0 = Vcjiimiic)) • ^°(m(Ca),m(6)) • 

where i = 0,... ,pa - 1, j = 0,... - 1, /c = 2,... ,pc, and (a, 6,c) = (1,2,3), (2,3,1), (3,1,2). There 

are a grand total of pq{r — 1) + qr{p — 1) + rp{q — 1) bubbles. 

6.4 Shape Functions 

As expected, they are the tensor products of the ID shape functions, and there are pqr such functions 
spanning 

6.4.1 Interior 

The coordinate free L? interior functions for the hexahedron are 

= [^’*](Mi(6))[^,](m(6))[^fc](m(6))(V/ii(ei)xVMi(6))-V;ui(6), (6.i7) 

for i = 0,..., p — 1, j = 0,..., (7 — 1 and /c = 0,..., r — 1. There are pqr interior functions. 

6.5 Orientations 

In 3D, both edges and faces have orientations associated to them, and they need to be considered to ensure 
full compatibility of shape functions along adjacent elements. Fortunately, these issues are handled almost 
effortlessly due to the structure of the formulas for the shape functions, and the local-to-global permutation 
functions, a^, and In what follows, it is assumed that §3.3, §4.5 and §5.5 have been covered. 

To construct orientation embedded shape functions, the first step is to predefine a sef of local orienfafions 
for each edge and face af fhe masfer elemenf level. After fhey are defined, fhese remain fixed. The nexf step is 
fo find fhe associafed locally ordered fuples of affine coordinafes represenfing fhose local orienfafions. Once 
fhese are found, fhe orienfafion embedded shape functions are merely fhe usual edge and face functions, buf 
wifh fheir respecfive ancillary operafor being precomposed wifh fhe appropriate local-fo-global permufafion 
funcfion. The only “burden” is fhen fo find fhe locally ordered fuples. This is shown nexf for fhe hexahedron. 

Figure 6.3 shows fhe masfer hexahedron along wifh a schematic represenfing all predefined local edge 
and face orienfafions. They represenf fhe o = 0 case. They are our choices for fhe local orienfafions (which 
in facf are fhe “lexicographic” orienfafions), buf ofhers may choose differenf local orienfafions fo represenf 
fheir o = 0 case. 
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Figure 6.3: Master hexahedron with numbered vertices and local edge and face orientations. 

To find the locally ordered tuples, the key is being aware of the relationships between the vertices and 
the affine coordinates. To illustrate this take as an example edge 12 and face 1234. 

Edge 12 is composed of the vertices ui and V 2 - Here, vi is linked to Ato(? 2 )^Ato(C 3 )^ while V 2 is 

linked to /Uo(^ 2 ) and /Uo(^ 3 ). The only difference between the the two vertices is that vi is linked to 

while V 2 is linked to |Ui(^i). Now, the local orientation is represented by the local vertex-ordering 
vi-^V2, SO quite simply the locally ordered pair is /roi(Ci) = (a*o(?i)j fio(?i)) (if the local ordering was 
V 2 -^vi, then the pair would be /rio(?i))- Hence, the orientation embedded edge 12 shape functions in 
with their gradient are 

= Fo(6)Fo(6)?^f (o-Q (/foi(6))), 

= f‘o(^c)/^o(6)V(/>f (fJo (/Ioi(6))) + </>f (o'o (/foi(6)))(/^o(6)V/io(6) + /^o(6)V/io(^3)) , 

where i = 2,... ,p. The same applies to the Tf(curl) edge 12 shape functions and their curl. Clearly the 
approach is analogous with any other edge. 

Face 1234 is composed of the vertices vi, V 2 , and U 4 . Here, the final goal is to find a locally ordered 
quadruple composed of two pairs. The local vertex-ordering corresponding to the local orientation of face 
1234 is vi-^V 2 -^V 3 -^V 4 . All one needs to do is to take the first two elements of the list, vi-^V 2 , and 
the second and third components of the list, namely V 2 -^vz- The former will represent the^zr^f pair in 
the quadruple, while the latter represent the second pair in the quadruple. Then one proceeds as if these 
where edges, so that vi-^V 2 is associated to /loi(li)! while V 2 -^vs is associated to / 2 oi(l 2 )- Finally the 
locally ordered quadruple is then the ordered succession of these two pairs, (/roi(Ci)) f(oi(C 2 ))- Hence, the 
orientation embedded face 1234 shape functions in with their gradient are 

= Aio(6)'^^(o-o (Foi(6),/foi(6))), 

= /^o(6)V</>°(cT°(/xoi(Ci),/ioi(6))) + </>*y(o-o (m(6),m(6)))V/io(?3) , 
where i = 2 ,... ,pa, j = 2 ,... ,pb, where {pa,Pb) are the orders of the first and second coordinate pairs in 
the quadruple cJq(/Ioi (Cl),/foi(? 2 ))- The same applies to the Ff (curl) andFF(div) face 1234 shape functions 
and their differential forms. Once again, the approach is analogous with any other face. 
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7 Tetrahedron 



Figure 7.1: Master tetrahedron with numbered vertices. 

The 3D simplex is the tetrahedron. The master element for tetrahedra is illustrated in Figure 7.1 in the 
X = {xi,X 2 ,xz) space. More precisely, it is the set {x G : xi > 0, X 2 > 0, X 3 > 0, xi + X 2 + X 3 < 1}. 

Denote vertex a by Va, so that vq = (0,0,0), vi = (1,0,0), V 2 = (0,1,0) and V 3 = (0,0,1). As 
described in §1.9, the 3D affine coordinates, Aq, Ai, A 2 and A 3 , are easily calculated for this master tetrahe¬ 
dron: 

Ao(x) = 1 - Xi - X 2 - X 3 , Ai(x) = Xi, X 2 (x) = X 2 , A 3 (x) = X 3 . (7.1) 

Their gradients are 

VAo(x) = (^- 1 ^ , VAi(x)=(^0, VA 2 (x)=(^i^, VA 3 (x)=(^ 0 . (7.2) 

These will be used explicitly or implicitly in what follows. 

Just like the triangle and the segment, the tetrahedron has a very natural correspondence of its vertices 
and its affine coordinates. Quite simply, each vertex Va is linked to the affine coordinate Aq, for o = 0,1, 2, 3. 
Indeed, Xa takes the value 1 at the associated vertex. 

Exact Sequence 

As with the hexahedron, the tetrahedron will have a 3D discrete polynomial exact sequence that represents 
the continuous exact sequence (1.7). It is 

pp ^ j\[P Tzq-P -pp-i ^ (7.3) 

where = 'PP(xi, X2, X3) is the space of polynomials of total order p. Meanwhile, the Nedelec and 
Raviart-Thomas spaces for the tetrahedron where already defined by (5.4) and (5.5), where N = Sin those 
definitions. 
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Like the triangle, the tetrahedron sequenee has an overall drop in polynomial order of one, which makes 
it compatible with the construction of the hexahedron. Also, as noted before, all of the spaces in the exact 
sequence are invariant under affine transformations. 

7.1 Shape Functions 

It will be clear that all the g(p + 3)(p + 2)(p + 1) shape functions lie in and span the space. 

The ideas in this section are completely parallel to those presented for the triangle (see §5) but in three 
dimensions. Therefore, the trace properties will not be analyzed in detail, since they follow analogously. 

7.1.1 Vertices 

The vertex shape functions and their gradients are simply the affine coordinafes fhemselves, 

(t)''[x) = \a{x ), V(j)'{x) = VAa(x), (7.4) 

for a = 0,1, 2, 3. There are a fofal of 4 verfex funcfions (one for each verfex). 


7,1.2 Edges 


These are freafed jusf like friangle edges. Hence, one can recur fo (l)f direcfly. Take for insfance edge 01. In 
Ibis case, fhe shape funcfions are simply 

^iix) = </>f(Aoi(x)) = (Ao(x) + Ai(x))>f 

'-^ --- 

'’‘“d project 

'-V- 

evaluate 


Ao(x) 


Ai(x) 


An(a;)+Ai(a:) ’ An(a;)+Ai(a:) 


for f = 2,..., p. The projection being implied is 


{xi,X2,X3) I )■ I )■ ) 0) 0) • 

If consisfs of finding fhe infersecfion P" = (, 0, 0) of fhe edge wifh fhe projecting plane passing 
fhrough fhe original poinf P = {xi,X 2 , x^) and fhe opposite nonadjacenf edge. If is illusfraled in Figure 7.2. 
Alfernafively if can be inferprefed in fwo sfeps. Firsf if is projected fo a poinf P' in an adjacenf face, using 
fhe projecfing line passing fhrough P and fhe disjoinf verfex fo fhe face. Once in fhe face, if is projecfed 
again fo fhe desired edge using fhe fradifional two dimensional friangle projecfion (see Figure 5.2). 

More generally, fhe edge functions and fheir gradienfs are 

<P^,ix) = 4>f{Aabix)), Vcj2‘r(x)=V4>f{Xab{x)), (7.5) 


wifh z = 2,...,p, 0<a<6<3. There are a fofal of p — 1 edge functions for every given edge, leading fo 
a fofal of 6(p — 1) edge functions. 
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Figure 7.2: Face projection from P to P' followed by an edge projection from P' to P". 

7.1.3 Faces 


The construction of these shape functions follows simply by homogenizing the triangle face bubbles, 
since this will represent a polynomial extension preserving the desired vanishing properties. This explains 
the definition of (j)fj in terms of homogenized polynomials. As an example, consider face 012, where the 
shape functions are 

= <^g(Aoi 2 (x)) = (Ao(x) + Ai(x) + A2(x))*+^;0g ( An(x)+Ai(x)+AAx) ^012(^) ) > 

-V- ^ 

project 

'-V-' 

evaluate 

for i > 2 and j > 1. The projection is already illustrated in Figure 7.2 and consists of finding fhe infersecfion 
P' = of fhe face wifh fhe projecfing line passing fhrough fhe original poinf P = {xi,X 2 ,x^) 

and fhe opposite verfex fo fhe face. 

The full collecfion of shape funclions and fheir gradienf is 

= 4>ij{^abc{x)) , = V(j)fj{Xabc{x)) , (7.6) 

where i > 2, j > 1, n = i + j = 3,... ,p, and 0<a<6<c<3. There are ^{p — l)(p — 2) shape 
funclions for each face, leading fo a fofal of 2 {p — l){p — 2) face funclions. 


7,1.4 Interior Bubbles 


The lelrahedron bubbles are given by blending a face shape function wifh a polynomial of complementing 
order which vanishes on fhe remaining face. As wifh Iriangles, if is carefully chosen as a Jacobi polynomial 


The inferior funclions and fheir gradienf are 

= 0g(Aoi2(*))[if+"’1 (moi(A3(i))) . 

V.#?,tW = |if+ti(m(A3W))V<Ag(X„i2M)+^g(A„2W)V[Lf+''](M0i(AsW)), 
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where i>2, j>\, k>l and n = i + j + /c = 4,. .. ,p, and where floi{\z{x)) = (1 — A 3 (x), A 3 (x)). 
There are g(p — l){p — 2)(p — 3) interior shape functions in total. 

7.2 i?(curl) Shape Functions 

The dimension of in three dimensions is ^p{p + 2){p + 3). A careful count of the linearly independent 
shape functions to be presented throughout this section will coincide with that dimension. Showing that the 
functions constructed are in follows from Lemma 2. The constructions are all analogous to those of the 
triangle and simply require of an extra extension which is naturally provided by homogenization. 


7.2,1 H (curl) Edges 


These are just the same as in the triangle case, but using three dimensional affine coordinates for the homog¬ 
enization. For example, for edge 01, the shape functions are 


Et{x) 


= E^{Xoi{x)) = [Pi](Aoi(x))(^Ao(x)VAi(x) - Ai(x)VAo(x 
= (Ao(x) + X^ix)ym 
= iXo{x) + Xiix)y+^Efl' 


Aifx) 


blend 


i \ \o(x)+Xi(x)^ Xo(x)+Xi(x) J ’ 

'■v'' 

project 


evaluate 


for f = 0,...,p — 1. Regarding the traces, note that they are completely inherited from F^^(Aoi(x)), which 
is a Whitney function known to have the desired vanishing properties and being tracewise compatible with 
the lower dimensional triangle edge functions. Therefore, all trace properties are satisfied, including fhe 
nonzero decay along fhe adjacenf faces fo fhe edge. 

The edge funclions wifh fheir curl are 


Efix) = EfiXabix)), V X Efix) = V X Ef (A,b(x)), 


(7.8) 


for i = 0,... ,p — 1, and 0 < a < b < 3. There are a fofal of p edge funclions for every given edge, for a 
lolal of 6p edge funclions. 


7,2,2 H (curl) Faces 

Like fhe Iriangle, fhe lefrahedron has Iwo families of shape funclions for every face. The Irace properlies 
follow from Ihose of Ihe edge funclions. There is a grand lolal of 4p{p — 1) face funclions. 
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Family I: The shape funetions and their eurls are 


Elj{x) = Ef-{Xabc{x)), V X Efj{x) = V X Efj{Xabc{x)), (7.9) 

for i > 0, j > 1, n = i + j = 1,..., p — 1, and 0<a<b<c<3. For every faee, there are ^p{p — 1) 
faee functions in this family. 

Family II: The shape functions and their curls are 

Ejjix) = Efj{Xbca{x)) , V X Ejj{x) = V X Efj{Xbca{x)) , (7.10) 

for i > 0, j > 1, n = i + j = 1,... ,p — 1, and 0<a<b<c<3. The only difference with the first 
family is that the entries are permuted to Xbca{x) instead of Xabc{x)- For every face, there are ^p{p — 1) face 
functions in this family. 

7.2.3 H (curl) Interior Bubbles 

The construction is completely analogous to that of in the sense that they are obtained by multiplying the 

2 f 2 ~1“ 7 ) 

face functions by the Jacobi polynomial Lf} . One must attempt this for various possible permutations of 
the entries, but being careful to ensure that they are linearly independent. Three families arise. 

The interior bubbles and their curl are 

VxE;=i.(x) = [Lf+''](/;„(Ai(a:)))VxEA(A.,,(i)) + V|Lf+'>](wi(Ai(i)))xEA(X„,„(a:)), 

where i > 0, j > 1, /c > 1, n = i+j + k = 2,... ,p— 1 and (a, b, c, d) = (0,1, 2, 3), (1, 2, 3, 0), (2, 3, 0,1), 
and where fiQi{Xd{x)) = (1 — Xd{x), Xd{x)). There is a grand total of \p{p — l)(p — 2) interior shape 
functions. 

7.3 i?(div) Shape Functions 

The dimension of TZT^ in three dimensions is ^p{p + 1) (p + 3). A careful count of the linearly independent 
shape functions presented here will coincide with that dimension. Showing that the functions constructed 
are in TZT^ is not immediate, but follows from the next lemma, which should be kept in mind. 

Lemma 3. Let x G M^/or N = 3, and fn G V^{x) be any polynomial of total order n in the coordinates 
X = {xi, X 2 , X 3 ). Given sq. si and S 2 affine coordinates in {or simply linear functions in x), it follows 
that the Raviart-Thomas space of order n + 1, contains the function 

fn{ • ) (^SqVsi X VS2 + SlVs2 X Vsq + S 2 VS 0 X Vsi^ G . 
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Proof. Recall the definition of the Raviart-Thomas space in three dimensions, 

TZT^ = 0 jy G (P^)^ :V{x) = f{x)x = with f G . 

Affine coordinafes are linear funclions in x = (xi, X 2 , 0 : 3 ), so fhaf 

Sfc(x) = Ok + bk-x, 

for ttk G M, G and k = 0,1, 2. Then Vsfc(x) = bk and 

V{x) = so(a^)Vsi(x) X 'Vs 2 {x) + si(x)Vs 2 (a;) x Vso(x) + S 2 (a^)Vso( 2 ;) x Vsi(x) 

= (oo + bo • x)( 6 i X 62 ) + (at + bi • x){b2 x 60 ) + (02 + • x)( 6 o x 61 ) 

= (^ao(bi X 62) + 01(^2 X 60) +a2(&o x 61)^ + (bo ■ {bi x b2)^x, 

'-V-" '-V-" 

=A =B(x) 

where fhe lasf ferm follows from various idenfifies. Clearly, A G (P°)3 = m 3 and bo ■ (bi x 62 ) G P° = M, 
so fhaf P G {P G : V(x) = f(x)x}. Hence, V G PT^. 

Now, fn € P” = P”“^ 0 P”, for n > 1 can always be decoupled info /„ = fn-i + fn, where 
fn-i G P”“^ and /„ G P”. As aresulf 

fn{x)V{x) = fn{x)A + fn-l{x)B{x) + Jn{x)B{x) , 

where if is clear /„A 0 fn-iB G (P"')3 and /„P G {V G (P”+i)3 ; l/(x) = (j){x)x}. Therefore, 
fnV G P7^+^ □ 

7.3.1 iT(div) Faces 

The general formula for fhese funclions is motivated by fhe well known firsl order Whilney form for H (div), 
along wilh fhe facf fhaf fhe normal frace of fhe faces should span fhe Iwo dimensional space. The general 
definilion is presenled nexf. 

Definition. Let sq. si and S 2 be arbitrary functions of some spatial variable in M^, with N = 3, and denote 
by p the order in the coordinate triplet (sq, S 2 ). Then 

Vif{so,si,S 2 ) = [P,-ff*^^](so,si,S 2 )(soVsi X VS 2 + SlVs 2 X Vso + S 2 VS 0 X Vsi^ , (7.12) 

fori = n—j, j = 0,... ,n and n = 0,... ,p—l {or equivalently f > 0, j > 0 and n = i+j = 0 ,... ,p—1 ). 
The divergence is 

V • si, S 2 ) = {i + j + 3)[Pi, Pj*'''^](so, 'Si 02 )Vso • (Vsi x VS 2 ) • (7.13) 

The formula for fhe divergence follows from fhe following lemma, because [Pi, Pj*^^](so, si, S 2 ) is a 
homogeneous polynomial of lolal order n = i + j in so, si and S 2 . 
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Lemma 4. Let ^n{so, si, ^ 2 ) £ 'P^iso, si, 'S 2 ) be a homogeneous polynomial of total order n in sq, si and 
S 2 , where so. si and S 2 are arbitrary functions of some spatial variable in with N = 3. Then 

V- (^V'n(so, Si, S2)(soVsi X Vs2+siVs2X VS0+S2VS0 X Vsi)^ = (n+ 3 )V’n(so, Si, S2)Vso-(Vsi X VS2) . 
Proof Let Lg'^(so, si, S2) = (sqVsi x Vs 2 + siVs2 x Vsq + S2VS0 x Vsi). First notice that 

V • (^so(Vsi X Vs2)^ = Vso • (Vsi X VS2) + so(^Vs2 • V X Vsi — Vsi • V x Vs2^ 

= Vso • (Vsi X VS2) , 

and similarly with V • (si(Vs2 x Vsq)) and V • (s2(Vsi x VS2)). All of them result in a scalar triple 
product, which is invariant to cyclic permutations. It follows 

^ ■ ^o^(' 5 o, Si, S2) = 3 Vso • (Vsi X VS2) . 

Now, consider a monomial SqsJs^. Then 

V(sgs5s2 ) • Loo(so,si,S 2) = (asg'^s^s^Vso + bsQs\~^S2Vsi +csgs5s2“^Vs2) • Loo (so, Si, S2) 

= OSq S^S 2 SoVso ‘ (Vsi X VS2) “t“ ^SqS^ S2S1VS1 ‘ (Vs 2 X Vsq) 

+ CSqs\s2~^S2Vs2 ■ (Vso X Vsi) 

= (ft b c)SqS]^S2Vso * (Vsi X VS 2 ) • 

With these last two results it follows 

V • (^SqS^s^Loo (so, si, 52)^ = (ft + 6 + c + 3 )sos 5 s 2 Vso • (Vsi x VS2) • 

Then observe that any homogeneous polynomial V'n(so, si, S2) is composed of monomials of the form 
SqsJs^ affixed total order a + b + c = n. The result immediately follows. □ 

The result in ( 7 . 13 ) is quite remarkable in the sense that, as required by the construction, there are 
no derivatives of [P*, P^^*'''^](so, si, S2) in the expression for the divergence. Although not used in this 
section, record the following useful remark which will be exploited when dealing with the prism and pyramid 
elements. 

Remark. Let vq = 1 — vi — V 2 , where ui and V 2 are arbitrary functions of some spatial variable in 
with N = 3, and where p is the order in the coordinates { 00 , 01 , 02 ). Then for all i > 0, j > 0 and 
n = i + j = 0,... ,p — 1, 

Vif{oo, oi, 02) = [Pi, Pf""^^]{oo, oi, 02 )Voi X V02 , V • Vif{oo, 01,02) = 0 . ( 7 . 14 ) 
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These functions are still represented by the same projection, and follow the logic of projecting, evaluat¬ 
ing and blending. For example, take face 012, so the shape functions are 


^ix) = ViiXoujx)) = jXoix) + Al(x) + A 2 (x))-+^'+^ Aoi 2 (x) ) , 


blend 


project 


evaluate 


for i > 0 and j > 0. It can be readily checked that they satisfy the required vanishing properties by focusing 
on the lowest order function 1 /qq(x). Moreover, by Lemma 3 it is clear VX G forn = i + j. 

More generally, the shape functions and their divergence are 


v!j{x) = VifiXaUx)) , V • (x) = V • V^fiXaUx)) , (7.15) 


where i > 0, j > 0, n = i + j = 0,... ,p — 1, and 0<a<6<c<3. Clearly, for each face there are 
^p{p + 1) functions, giving a total of 2 p{p + 1) face functions. 


7.3.2 H (div) Interior Bubbles 

The construction is just like in and Tf(curl), and there are three resulting families to consider. The 

2 f 2 7 ~1“ 1 j 

polynomial chosen to vanish is the Jacobi polynomial, 

The interior functions with their divergence are 

where i > 0, j > 0, /c > 1, n = i+y + A; = 1,... ,p— 1 and (a, b, c, d) = (0,1, 2, 3), (1, 2, 3, 0), (2, 3, 0,1), 
and where jxoi{Xd{x)) = (1 — Xd{x), Xd(x)). There is a grand total of ^p(p — l)(p + 1) interior shape 
functions. 


7.4 Shape Functions 

The gp(p + 2)(p + 1) shape functions will clearly span PP~^. 

7.4.1 Interior 

Similar to the construction of triangle functions, the interior functions for the tetrahedron are 

= [Pi|(A„(a,).Ai(i))|F,!“+‘](A„(j,) + Ai(j,).A2(a,))|pf «+''|(1-A3W,A3 (,t)) 

= [Pi,P/‘+'l(Aoi2(i))|Pf+"+"1 (wi(A3(i)))(VAi{x)xVA2W)-VA3(i) , 

where i > 0, y > 0, A: > 0, and n = i + j + k = 0,...,p —1, and where poiiXsix)) = (1 — A 3 (x), A 3 (x)). 
There is a total of gp(p + 2)(p + 1) interior functions. 
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7.5 


Orientations 



Figure 7.3: Master tetrahedron with numbered vertices and local edge and face orientations. 

To construct orientation embedded shape functions for the tetrahedron, it is recommended to have read at 
least the first part of §6.5. The predefined local edge and face orientations for the tetrahedron are illustrated 
in Figure 7.3. They represent the o = 0 case. The task at hand is to find the associated locally ordered tuples 
of affine coordinates representing those local orientations. As examples take edge 01 and face 012. 

Edge 01 is composed of the vertices vq and vi, which are linked uniquely to Aq and Ai respectively. 
The local orientation for edge 01 is represented by the local vertex-ordering no ni. Asa result, the locally 
ordered pair for edge 01 is Aqi = (Aq, Ai). It follows that the orientation embedded edge 01 shape functions 
m with their gradient are 

= (/>f ((t^(Aoi(x))) , V(/>|(x) = V(/)f (fj^(Aoi(x))), 

with i = 2,... ,p. The same applies to the H (curl) edge 01 shape functions and their curl. The approach is 
analogous with any other edge. 

Face 012 is composed of the vertices no, ni and V 2 , which are linked to Ao, Ai and A 2 respectively. The 
local orientation for the face is represented by the local vertex-ordering vo-^vi-^V 2 , and as a result the 
locally ordered triplet for face 012 is A 012 = (Ao, Ai, A 2 ). Hence, for example the orientation embedded 
face 012 shape functions with their gradient are 

= (/)g(a^(Aoi2(x))), V0^j(x) = V</>g(o-f^(Aoi2(x))), 

with i > 2, i > 1 and n = i + j = 3,... ,p. The same applies to the F7(curl) and i7(div) face 012 shape 
functions and their differential forms. Naturally, the approach is analogous with any other face. 
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8 Prism 




Figure 8.1: Master prism with numbered vertices. 

The master prism is shown in Figure 8.1 in the (x, z) = (xi, X2, z) space. It is the Cartesian product of 
a triangle and a segment. More specifically, it is {x G : xi > 0, X2 > 0, xi + X2 < 1} x (0,1). Here, x 
represents the triangle element coordinates, while z represents the ID segment element coordinate. 

The prism has both quadrilateral and triangle faces. Similarly, there are two types of edges. Those edges 
which are adjacent to both a quadrilateral face and a triangle face are called mixed edges, while those edges 
only shared by quadrilateral faces are simply called quadrilateral edges. These distinctions are important, 
and the form of the shape functions will differ for the different types of edges and faces. 

Due to the Cartesian product structure, it is natural to consider the 2D affine coordinates for the triangle 
(dependent on x) and the ID affine coordinates for the segment (dependent on z). For this master element 
they are 

^^{x) = 1 — Xl — X 
do{z) = 

Their gradients in 3D are 

Vuo{x) = (-i) > 

Vtioiz) = 

These are important and are used explicitly or implicitly in many of the oncoming calculations. 

As usual, there are natural relationships between vertices, edges and faces, and the affine coordinates. 
The related coordinates are those which take the value 1 at the given topological entity. In fact, for the 
prism each vertex is linked to two affine coordinates, one of them a 2D affine coordinate and the other a ID 
affine coordinate. Moreover, each edge is linked to one affine coordinate. For mixed edges it is a ID affine 
coordinate, while for quadrilateral edges it is a 2D affine coordinate. Lastly, the triangle faces are also linked 
to one ID affine coordinate. For example, vertex 0, vq = (0,0,0), is linked to the affine coordinates oo{x) 


2, Ui{x)=Xi, I/2(x)=X2, 
1- Z, Hl{z) = z . 


(8.1) 


Vz^i(x) = 


Vl'2{x) = 


V/ri(2;) = 0 


( 8 . 2 ) 
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and fio{z). Meanwhile, mixed edge 01 is linked to the affine eoordinate fio{z), while quadrilateral edge 03 
is linked to the affine eoordinafe Finally, faee 012 is linked fo ixq{z). 


Exact Sequence 

The produel slruelure of Ihe prism suggesls lhaf Ihe diserele polynomial exael sequenee approximaling (1.7) 
is inlimalely relaled lo Ihe diserele sequenees for Ihe Iriangle and segmenl (see (5.3) and (3.3)). Indeed, Ihis 
is Ihe ease. The sequenee is of Ihe form 

Vpp:-? QP.? yp.9 ^ (8.3) 

where 

= VP = VP{xi,X2)®V‘i{z), 

QP’^ = (AfP ® p9) X {VP ® V^-^) , 

(8.4) 

= (7^rP ® p^-i) X {VP-^ ® V^) , 

yp,g ^ -pp-l ^ pq-l ^ 

Here, Ihe spaees AfP and TZTP’ eorrespond lo Ihe N = 2 ease (see (5.4) and (5.5)), so lhal Ihe spaee 
AfP (g) 'pt = {(j)E : (f) G V^, E G AfP} has Iwo eomponenls and dimension p{p + 2)(q + 1). The same 
applies lo TZTP <8) V^~^ whieh has dimension p{p + 2 )q. 

8.1 Shape Functions 

It will be elear that all shape funetions lie in VP®V^, whieh has dimension i(p+2)(p+l)(q+l). Moreover, 
a judieious eount of the shape funetions eonstrueted in this seetion eoineides with that dimension, ensuring 
that the span of these is preeisely (g P'^. 

Notiee that for the triangle there are vertex, edge and faee shape funetions, while for the segment 
there are vertex and edge funetions. The six possible tensor produets of these will preeisely give the 
shape funetions for the prism. The vanishing properties are naturally inherited from eaeh of the eomponents 
of the tensor produet strueture of the shape funetions, so they follow easily. 

8.1.1 Vertices 


As usual with Cartesian produet struetures, the vertex funetions are simply the produet of the affine eoordi- 
nates assoeiated to the given vertex. Henee, they are the tensor produet of lower dimensional vertex funetions 
whieh inherit all the desired vanishing properties and the deeays along adjaeent faees to the vertex. 

The vertex shape funetions and their gradients are 

(IA{x,z) = Va{x)pc{z), 

(^O. J J 

V(IA{x, z) = Pc{z)Vl'a{x) + l'a{x)Vpc{z) , 
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where a = 0,1, 2 and c = 0,1. Notice c = 0 represents the “bottom” face (vertices 0, 1 and 2), while c = 1 
represents the “top” face (vertices 3, 4 and 5). Also Va represents vertices a and a + 3. There is a total of 6 
vertex functions (one for each vertex). 


8.1.2 Edges 


Mixed Edges. These functions are the tensor product of triangle edge functions with ID 
vertex shape functions. For instance, take the edge 01, which is in the bottom triangle face (associated to 
/io(^))- The shape functions then take the form 


(fiix) = lXQ{z)4)f{vQl{x)) = hq{z){vo{x) + Vl{x)Y (t)\ 






blend 


Uo(x)+I/l(x}'> l/o(x) + l'l(x) 

' -V-' 

project 


evaluate 


for i = 2,..., p. The trace properties are naturally inherited along the edge, its adjacent faces (including the 
nonlinear decay in the triangle face when = 1) and all the other faces where it is supposed to vanish. 
Like the triangle, the projection is of the form 


{XI,X2,Z) 


{XI,X2,0) 




It consists of finding the intersection P" = 0,0) of the edge with the plane passing through the orig¬ 

inal point P = {xi,X 2 , xs) and the opposite disjoint quadrilateral edge. Alternatively, it can be interpreted 
as a two step projection, where it is first projected to an adjacent face, and then projected to the desired edge. 
This projection is shown in Figures 8.2 and 8.3. 




Figure 8.2: Triangle face projection from P to P' followed by a mixed edge projection from P' to P". 


The shape functions with their gradient are 

())|(X, Z) = ^c{z)(l)f{Uab{x)) , 

V())|(X, Z) = + 4)f{iJab{x))V^ic{z) , 
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where i = 2,...,p, 0<a<b<2, and c = 0,1. There are p — 1 shape functions for each mixed edge, for 
a total of 6 (p — 1 ) mixed edge functions. 


Quadrilateral Edges. These are the tensor product of triangle vertex shape functions with the ID 
edge functions. For edge 03, the shape functions are 

(^|(x, z) = yo{x) c/)f{ po{z),^pi{z) ), 

blend project 

^-V-^ 

evaluate 

for i = 2,... ,p. As expected, there is a linear edge blending towards both of the adjacent quadrilateral 
faces, given by while all other trace properties are also inherited. The implied projection is simply 


{xi,X2,z) I—^ 


i 0 , 0 ,z). 


It consists of finding the intersection P'" = (0, 0, z) of the edge with the normal plane passing through 
the original point P = (xi, X 2 , X 3 ). Altemaltively it can be interpreted as a two step projection. This is 
illustrated in Figure 8.3. 

i i 


Quadrilateral Face Projection Edge Projections 

Figure 8.3; Quadrilateral face projection from P to P' followed by a mixed edge projection from P' to P" and a 
quadrilateral edge projection from P' to P'". 

The shape functions and their gradient are 




z) = Va{x)(l)f{llm{z)) , 

V(?!)|(x, z) = Va{x)V(l>f{flm{z)) + (t)f{floi{z))Vva{x ), 


where i = 2,... ,q and a = 0,1, 2. There are q — 1 shape functions for each quadrilateral edge, for a total 
of 3{q — 1) quadrilateral edge functions. 
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8.1.3 Faces 


Triangle Faces. These are tensor products of triangle face bubbles and ID vertex shape 
functions. For instance, for face 012 the shape functions are 

z) = //o(^) 1^1 (x), 1 / 2 (x)) , 


blend 


project 


evaluate 

for i > 2 and j > 1. The trace properties are trivially inherited from each of the components. The triangle 
face projection is illustrated in Figure 8.2 and consists of finding the intersection P' = (xi,X 2 , 0 ) of the 
face with the normal line passing through the original point P = {xi,X 2 , X 3 ). 

In general, the shape functions and their gradient are 

4>\j{x,z) = 


( 8 . 8 ) 


'^(Pijix,z) = ^c(^)V(/>j)(Foi2(x)) + (pij{i7oi2ix))V , 

where i > 2, j > 1, n = i + j = 3,... ,p, and c = 0,1. As with any triangle face, there are 
^{p — l)(p — 2) shape functions per face, for a total of {p — l){p — 2) triangle face functions. 


Quadrilateral Faces. The quadrilateral face functions are simply the product triangle edge func¬ 
tions and ID edge functions, so they conveniently fall into the general definition of . For face 0143, 
they take the form 


4’lj{x,z) = <p°j{Poi{x),fIoi{z)) = (;/o(x) + z/i(x))>g( poijz) 


blend 


project 


evaluate 

where i = 2,... ,p and j = 2,... ,q. Clearly, the desired trace properties are inherited. The implied 
projection, already illustrated in Figure 8.3, consists of finding fhe infersecfion P' = 0, z) of fhe face 

wifh fhe projecting line lying in fhe horizonfal plane and passing fhrough fhe original poinf P = (xi, X 2 , X 3 ) 
and fhe opposite disjoin! quadrilateral edge. 

The shape functions and fheir gradienf are 

(t^ij{x,z) = 4>fj{iyab{x),j2oi{z )), V(/)-y(x, 2 ;) = V(/>° (f4fe(x), ^oi(^)), (8.9) 

where i = 2,... ,p, j = 2,... ,q, and 0 < a < b < 2. There are {p — l){q — 1) shape functions per face, 
leading fo a fofal of 3{p — l){q — 1) quadrilateral face functions. 


8.1.4 Interior Bubbles 

The bubble functions are the tensor product of the triangle face functions and ID edge functions. 
Due to their structure, the zero trace properties are trivially satisfied along the whole boundary. 
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The interior bubbles and their gradient are 

'^4>^jkix,z) = (j)kim{z))V<j)f^j{uoi2{x)) + (j)^j{uoi2{x))V4>^{poi{z)) , 

where i > 2, j > 1, n = i + j = 3,... ,p, and k = 2,... ,q. There are a total of ^{p — l)(p — 2){q — 1) 
bubble shape funetions in H^. 

8.2 i?(curl) Shape Functions 

The linearly independent shape functions presented here are shown to belong and span the H (curl) con¬ 
forming space {AfP (g) V^) x {V^ ® which has dimension p(p + 2)(g + 1) + |(p + 2)(p + l)q. 

These shape functions, as expected, are composed of combinations of and Tf(curl) components. 
Intuitively, they involve the affine coordinates and at least Ef, Efj and E^. All shape functions continue to 
respect the logic of projecting, evaluating and blending, and for a given topological entity, the projections 
are the same as those in (see Figures 8.2 and 8.3). 


8.2,1 H (curl) Edges 


H (curl) Mixed Edges. These are tensor products of triangle H (curl) edge functions and ID vertex 
functions. For example, take edge 01. Note that Ef{vQi{x)) in three dimensions is a three component 
vector whose last component is zero, since it is independent of the z coordinate. Indeed, if considered 
in two dimensions, it is just the well defined friangle Ff(curl) edge funcfion. As such, if is an elemenf 
of NP (in fwo dimensions). Due fo fhis imporfanf facf, fhe firsf fwo componenfs of Ef are denofed by 
{Ef{vQi{x)))x G For fhis edge, fhe shape funcfions are 


^ PQ{z){Ef{vm{x)))^ 
0 


Et{x,z) = pi,{z)Ef{Voi{x))= 




blend 


project 


evaluate 


for i = 0,... ,p — 1. Hence, E® G {AfP (g) 'P'^) x {VP (g) V^~^). For fhe frace properfies, nofe fhaf over faces 
1254 and 0253 bofh fangenfial componenfs are zero, since fhe z componenf is zero and fhe fangenf fo edges 
12 and 02 is zero by inherifance of Ef{vQi{x)). In fhe lop face, fhe fangenf is also zero because pq{z) is 
zero fhere. Finally, by consfrucfion, over face 012 ifs fangenf is precisely fhe friangle iF(curl) edge shape 
funcfion {Ef (Fqi {x)))x, while over face 0143 ifs fangenf behaves like a quadrilateral H (curl) edge function 
l^o{z){Ef{fIoi{xi)))x. 
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Now, the general shape funetions and their eurl are 


Et{x,z) = lic{z)Ef{Vab{x)) , 

V X El{x, z) = Hc{z)V X Ef{Uab{x)) + V^c{z) X EfiUabix)) , 

with r = 0,...,p — 1, 0<a<6<2 and c = 0,1. There are p shape funetions for eaeh mixed edge, giving 
a total of 6p mixed edge funetions. 


iT(curl) Quadrilateral Edges. These are tensor produets of triangle vertex shape funetions and 
Ef{poi{z)), whieh aet as ID i7(curl) funetions (even though these do not formally exist for the segment 
element). For instanee, take edge 03. Now, Ef{ilQi{z)) in three dimensions is a three eomponent veetor 
whose first two eomponents are zero, sinee it is only dependent on the z eoordinate. Indeed, if eonsidered 
in one dimension, it is just a segment funetion belonging to This last eomponent of Ef is denoted 

by {Ef{flQi{z)))z G V^~^. For this edge, the shape funetions are 


Et{x,z) 


blend project 


evaluate 


0 G {0} X {0} X (pi ® , 

\uQ{x){Ef{pm{z)))z) 


for f = 0,..., g — 1. Clearly Ef G ® V^) x (P^ (g) P'^“i). The traee properties follow easily for the 
triangle faees, sinee the veetor field points normal to those faees, while in the nonadjaeent quadrilateral faee 
1254 it is zero due to 1 ^ 0 ( 3 ;) being zero. In the adjaeent faees, Ef{jloi{z)) is unaffeeted by the restietions (it 
is already tangent to the faees and independent of x), so the tangential traee is a quadrilateral H (eurl) edge 
funetion simply beeause the restrietion of vq{x) is a linear blending funetion over that faee. 

In view of (4.10), the shape funetions and their eurl are 


Ef{x,z) = Ua{x)Ef{fIoi{z)), 


V X Ef{x, z) = Vva{x) X Ef{Jloi{z )), (8.12) 


with i = 0,..., g — 1, and a = 0,1, 2. There are g shape funetions for eaeh quadrilateral edge, for a total of 
3g quadrilateral edge funetions. 


8.2.2 H (curl) Faces 

H (curl) Triangle Faces. These are tensor produets of triangle H (eurl) faee funetions and vertex 
shape funetions. As usual, there are two families. Proeeeding as with H (eurl) mixed edge funetions, it 
follows the shape funetions lie in {J\f^ ® V^) x (P^’ ® and that the traee properties are satisfied. 

As aeeusfomed, there are p{p — 1) funetions per triangle faee, and a grand total of 2p{p — 1) triangle faee 
funetions. 
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Family I: The shape funetions and their eurl are 


Elj{x,z) = ^c{z)Ef-{voi 2 {x )), 

V X eI^{x,z) = lJLc{z)V X Ef-{vQi2{x)) +Viic{z) X Ef-{vQi2{x)) , 

for i > 0, j > 1, n = i + j = 1,... ,p — 1, and c = 0,1. For every faee, there are \p{p 
this family. 

Family II: The shape funetions and their curl are 

Elj{x,z) = pc{z)Ef-{Vi2o{x)) , 

V X eI^{x, z) = Pc{z)v X E^^{vi 2 o{x)) + Vpc{z) X E^j{vi 2 o{x )), 

for i > 0, j > 1, n = i + j = 1,... ,p — 1, and c = 0,1. Note the fact that the entries are Fi 2 o(a;) as 
opposed to Foi 2 (x). For every face, there are ^p{p — 1) functions in this family. 

H (curl) Quadrilateral Faces. These are tensor products of H (curl) triangle edge functions and ID 
edge shape functions, and viceversa, so that it is clear there are two families. Both fall naturally into the 
general definition of E^. There are a total of p{q — 1) + q{p — 1) functions per quadrilateral face, for a 
grand total of 3{p{q — 1) + q{p — 1)) quadrilateral face functions. 

Family I: To begin, take for example face 0143. The shape functions for the first family are of the 
form £'R(Foi(x), j2oi{z)), so proceeding exactly as with the mixed edges it is easily shown that it lies in 
X {0} and more importantly that the trace properties hold. 

The general shape functions and their curl are 

Elj{x,z) = Efj{vab{x),jlm{z)), V x £'^(x, z) = V x Efj{vab{x),floi{z)), (8.15) 

for i = 0,... ,p — 1, J = 2,..., q, and 0 < a < 6 < 2. There are p{q — 1) shape functions in this family 

per quadrilateral face. 

Family II: Again, taking face 0143 as an example, the second family of shape functions has the 
form E^{p,Qi{z), Fbi(x)). This time, proceeding as with quadrilateral edges, one can show that the trace 
properties hold and that the functions lie in {0} x {0} x {V^ (8) 

The general shape functions and their curl are 

Elj{x, z) = E^j{fIoi{z),Uab{x)), V xEIj{x,z) = V X E^j{floi{z),i7ab{x)), (8.16) 

for i = 0,..., q — 1, j = 2,... ,p, and 0 < a < 6 < 2. Note the entries are permuted with respect to the 

first family. There are q{p — 1) shape functions in this family per quadrilateral face. 


(8.13) 
— 1) functions in 
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8.2.3 H (curl) Interior Bubbles 


There are three families of interior bubble funetions. They involve tensor produets of 4>^, Ef and 
Efj. The first two families have elements lying in {AfP (g) 'P'?) x {0}, while the last family has elements in 
{0} X {0} X (PP (g) P”? ^). The traee properties are also satisfied by using similar argumenfs fo fhose used 
for fhe faee and edge funetions. There is a grand fofal of p{p — 1) (g — 1) + ^ (p — 1) (p — 2)g inferior bubble 
funetions. 


Family I: The shape funetions and fheir eurl are 

Eijkix, z) = (poi(^))^*y (z^ 0 i 2 (a:)), 

V X Efjk{x,z) = (l)^{poi{z))V X P^(z7oi2(x)) + V(l)^{f2oi{z)) x E^{Pou{x)), 

wifh i > 0, j > 1, n = f + j = 1,... ,p — 1, and k = 2,... ,q. The family has ^p{p — l){q — 1) funetions. 

Family II: The shape funefions and fheir eurl are 

Efjkix, z) = (l)l{poi{z))E^j{vi 2 o{x )), 

V X E^^^{x,z) = 4)l{poi{z))V X P^(i7i2o(x)) + V(j)l{iloi{z)) X P,^(z7i2o(a:)), 

wifh i > 0, j > 1, n = i + j = 1,... ,p — 1, and k = 2,... ,q. The only differenee wifh fhe firsf family is 
fhaf fhe permufafion iyuoix) is used insfead of iyoi 2 {x)- The family has ^p{p — l){q — 1) funefions. 


Family III: Due fo (4.10), fhe shape funefions and fheir eurl are 

Efjk{x,z) = (/>g(l7oi2(x))Pf(poi(^)) , 
u A t:. (8-19) 

V X Efj^{x,z) = V0g(i7oi2(x)) X Ef{iloi{z)), 

wifh i > 2, j > 1, n = f + j = 3,... ,p, and k = 0,... ,q — 1. The family has ^(p — l)(p — 2)g funefions. 

8.3 i?(div) Shape Functions 

The linearly independent shape funetions presented here are shown to belong and span the H (div) eonform- 
ing spaee {TZT^ (g) P*?"^) x (pP“^ (g) P*?), whieh has dimension p(p + 2)q + \p{p + l)(g + 1). 

As expeeted, the shape funetions are eomposed of eombinations of the lower dimensional H^, P(curl) 
and H (div) eomponents. Intuitively, they involve the affine eoordinates and at least and V^. Again, the 
shape funetions respeet the logie of projeeting, evaluating and blending, and for a given faee, the projeetions 
are the same as those in (see Figures 8.2 and 8.3). 
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8,3,1 H{dW) Faces 

H (div) Triangle Faces, These are the tensor product of H (div) triangle face functions and ID vertex 
functions. The Tf(div) triangle functions are of the form Vj^{i7oi2{x)), so hy (7.14), it follows that their 
first two components are zero. The last component, which corresponds to the normal trace, by construction 
is an triangle face function [Pi, P^]{vqi 2 {x)), which is known to lie in 'PP~^. For face 012, the shape 
functions have the form ii,Q{z)V^^{vQi 2 {x)), meaning that they lie in {0} x {0} x (pp~^ (g) V^). The trace 
properties also follow since the function points normal to the triangle faces, meaning that they are tangent 
to the quadrilateral faces, an as a result have zero normal trace along those faces. Meanwhile at the opposite 
triangle face 345, the function is also zero because iio{z) is zero there. Lastly, at the face itself, iio{z) is 
unity, so along the face the normal component is precisely an L? triangle face function. 

Due to (7.14), the triangle shape functions and their divergence are 

yljix, z) = ^c{z)Vif{voi 2 {x)), V • 1/^(x, z) = Vuciz) ■ V;^(Foi 2 (x)) , (8.20) 

where i > 0, j > 0, n = i + j = 0,... ,p — 1 and c = 0,1. There are ^p{p+ 1) shape functions per triangle 
face, leading to a total of p{p + 1) triangle face functions. 


H{dW) Quadrilateral Faces, These are cross products of 77(curl) edge functions. They naturally fall 
into the definition of V^. For instance, take face 0143. Now, FJf (Foi(x)) in 3D is a three component vector 
whose last component is zero, since it is independent of the z coordinate. Similarly, Ef{ilQi{z)) in 3D is a 
three component vector whose first two components are zero, since it is only dependent on the z coordinate. 
It then makes sense to speak of {Ef{i7Qi{x)))x G A/"^ and [Ef {flQi{z)))z G Moreover, as discussed 

in §5.3, ( Q )(£'j^(i7oi(x)))a; G For this face, the shape functions can be written as 


yij{x,z) = V^{voi{x),ilQi{z)) = 




(£;f(Foi(x))), 


0 


(i7f(^oi(^))). 


{Ef{n))z[M]{Ef{P^^)) 


G (7^7^®iP‘?■^) X {0} 


for f > 0 and j > 0. Hence, Vlj{x,z) G [TZT^ (g) x [pP~^ (g) 'P'?). The trace properties are 

also satisfied because Ef{j2Qi{z)) is normal fo fhe friangle faces, so 47j(x, z) musf be fangenfial, and as a 
resulf has zero normal frace. Af fhe quadrilaferal faces one only has fo look af fhe fangenfial componenf of 
Ef{^oi{x)) along fhe mixed edges, and if follows fhaf fhe normal fraces of Vlj{x, z) fo fhe faces 1254 and 
0253 are zero, while af face 0143 if fakes fhe form of an quadrilaferal face function. 

The quadrilaferal face functions and fheir divergence are 


yij{x,z) = V^{vab{x),iloi{z)), V • V^j{x,z) = V • V^{i7abix),floi{z)), (8.21) 


where i = 0,... ,p — l, j = 0,... ,q — l, and 0 < a < 6 < 2. There are pq shape funclions per quadrilateral 
face, for a fofal of 3pq quadrilaferal face funclions. 
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8.3.2 H (div) Interior Bubbles 


There are three families of -ff(div) bubble funetions. Two of them are closely related to and have 
elements in {TZT^ (g) 'P'J-i) x {0}, while the third family is related to an has elements in the space 
{0} X {0} X ^ (g) 'P'?). Naturally, the trace properties are satisfied by using similar arguments to those 
used for the face functions. There is a grand total of p{p — l)g + ^p{p + 1) (g — 1) interior bubble functions. 

Family I: Using (4.10), the shape functions and their divergence are 

Vijkix,z) = X Ef{floi{z)), 

u ^ . \ (8-22) 

V • = E^imiz)) • (V X , 

where i > 0, y > 1, n = i + y = 1,... ,p — 1, and A: = 0,..., g — 1. The family has ^p{p — l)q functions. 


Family II: Using (4.10), the shape functions and their divergence are 

V^jkix,z) = Ef^{i7i2o{x)) X Ef{poi{z)), 

. ^ . \ (8-23) 

V • Ub,(x,z) = E^ifloiiz)) • (V X , 

where i > 0, j > 1, n = i + j = 1,... ,p — I, and A: = 0,..., g — 1. The only difference with the first 
family is that the permutation 1 ^ 120 ( 2 ;) is used instead of uqi 2 {x)- The family has ^p{p — l)q functions. 


Family III: Due to (7.14), the shape functions and their divergence are 

y^jkix,z) = (j)k{miz))V^{izoi2ix)), 

V • z) = V<J)^{iloiiz)) ■ V^ii^ouix)), 

with i > 0, j > 0, n = f + j = 0,... ,p — 1, and k = 2,... ,q. The family has ^p{p + 1)((7 — 1) functions. 


(8.24) 


8.4 Shape Functions 

In this case, the space (g) P'S is spanned by the ^p{p + !)<? linearly independent shape functions. 

8.4.1 Interior 

These are the tensor products of L? triangle functions and edge functions. They are 

i^^jkix^z) = [Pi]{uo{x),vi{x))[Pf"^\uo{x) + ui{x),V2{x))[Pk]{Mz)-,lJ‘i{z)) 

(8.25) 

= [P,Pf+i](z7oi2(x))[Pfc](/7oi(^))(Vivi(x)xVi.2(x))-V/ii(z), 

where f > 0, j > 0, n = i + y = 0,...,p — 1 and A: = 0,..., g — 1. There is a total of \p{p + 1)(7 interior 
functions. 
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8.5 Orientations 



Figure 8.4; Master prism with numbered vertices and local edge and face orientations. 

To construct orientation embedded shape functions for the prism, it is recommended to have read §6.5 
and §7.5. The predefined local edge and face orientations for the prism are illustrated in Figure 8.4. They 
represent the o = 0 case. The task at hand is to find the associated locally ordered tuples of affine coordinates 
representing those local orientations. As usual, the key is being aware of the relationships between the 
vertices and the affine coordinafes. As examples lake edges 01 and 03, and faces 012 and 0143. 

For mixed edge 01, fhe vertices are vq and vi, which are linked to {vq{x) , ^q{z)) and {vi{x) , iJiQ{z)) 
respectively. The only difference is that vq is linked to t'o(x), while is linked to vi{x). The local 
orientation for edge 01 is represented by the local vertex-ordering uo-->^ui. Therefore, the locally ordered 
pair for edge 01 is Foi(x) = {vq{x), oi{x)). The orientation embedded shape functions for mixed edge 01 
are simply the usual shape functions, but with their respective ancillary operator and differential form (that 
is (j)f, Ef, V(/)f and V x Ef) being precomposed with and evaluated at the locally ordered pair. 

For quadrilalteral edge 03, composed of vertices vq and vs, the differing affine coordinates are fj,o{z) 
and (z) respectively. Since the local vertex-ordering is uq --)• us, it follows the locally ordered pair for this 
edge is p,Qi{z). Then, the orientation embedded shape functions are constructed like those of mixed edge 
01. That is, precomposing the ancillary operators with and evaluating at the locally ordered pair. 

For triangle face 012, composed of vertices vq, and V2, the differing affine coordinates are t'o(x), 
r'i{x) and i'2{x) respectively. The local vertex-ordering is vo-^vi-^V2, so the locally ordered triplet for 
this edge is Foi 2 (x). Then, the orientation embedded shape functions are the usual shape functions but with 
the ancillary operators Efj, Vj^, V x Efj and V • Vj^) precomposed with and evaluated at 

the locally oriented triplet. 

Finally, quadrilateral face 0143 has local vertex-ordering vq-^vi-^v^-^vs, so one only needs to look 
at uo--)'Ui and vi-^v^ as if they were edges. This leads to the locally ordered pairs Foi(x) and iloi{z) 
respectively, so the locally ordered quadruple is (Foi(x), floi{z)). Again, the orientation embedded shape 
functions are simply the shape functions but with the ancillary operators ((^?-, E^, V^, V(/>?-, V x 77^ and 
V • V^) precomposed with cj^ and evaluated at the locally oriented quadruple. 
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9 Pyramid 



Figure 9.1: Master pyramid with numbered vertices. 

The master pyramid is shown in Figure 9.1 in the (.^,C) = space. More specifically, the 

definition is {(6,6,0 € : 6 > 0 ,6 > 0, C > 0 ,6 + C < 1 ,6 + C < 1}- 

Clearly, the pyramid is neither a simplex nor a Cartesian product, but it captures features of both quadri¬ 
laterals and triangles. Indeed, it has both quadrilateral and triangle faces. Similarly, it has two types of 
edges. Those edges which are adjacent to both a quadrilateral face and a triangle face are called mixed 
edges, while those edges only shared by triangle faces alone are called triangle edges. These distinctions 
are fundamental, and the form of the shape functions will differ for the different types of edges and faces. 

Due to the virtually unknown structure of the pyramid, at first it seems almost like an insurmountable task 
to be able to find representative functions that resemble affine coordinates for this 3D element. Surprisingly, 
there is in fact such a set (or rather set^) of coordinates. However, to reach that point, it is better to start 
by elementary means. With this in mind, the idea is to separately analyze the affine coordinafes of the 
quadrilateral face and the triangle faces. 

In truth, nothing is getting in the way of explicitly computing the 2D triangle affine coordinates of each 
of the four triangle faces as described in §1.9. There turns out to be two independent sets of such coordinates 
which are 

i^o( 6 ,C) = 1 -6 -C, z^i( 6 ,C) = 6 , 1 ^ 2(0 = C, 

1 ^ 0 ( 6 ,C) = 1 -6 -C, z^i( 6 ,C) = 6 , 1 ^ 2(0 = C- 

Indeed, the triplet Foi 2 ( 6 , C) represents triangle faces 125 and 435, while the triplet ^ 012 ( 6 , C) represents 
triangle faces 145 and 235. Their gradient is 

Vz^o(6,C) = , Vz/i(6,C) = ( 1 ), 

Now, the quadrilateral face can undergo a similar treatment, resulting in the standard two sets of ID 


Vz^o(6 ,C)= 0 


Vz^i(6,C) = 


Vi^2(C) = 
Vi^2(C) = 


(9.2) 
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affine coordinates, 


/^o(Ci) = i-6> f^i(6) = 6, 

/^o(?2) = i- 6> f^i(6) = 6- 

However, these are convenient only when restricted to the 2D quadrilateral face, and not in 3D. The reason 
is that they do not act as blending functions to the. faces. For example, /ro(^ 2 ) is unity at face 125, but it does 
not vanish at the opposite face, which is face 435. This inconvenience does not occur with the hexahedron 
or prism due to the Cartesian product structure of those elements. Despite this setback, it is possible to fix 
this by considering scaled coordinates which additionally depend on C,. The sets of quadrilateral scaled ID 
affine coordinates are 

= 1 - 1^ , = 

These can be readily checked to act as face blending functions between the opposite triangular faces, which 
is precisely what was desired. Moreover, when restricted to the quadrilateral face, so ^ = 0, they coincide 
with the usual sets of affine coordinates for the 2D quadrilateral faces. Their gradient is 


VFo(i^)= (t4^( 

Vfo(i^) = ^(-(W)), 




(9.4) 


Lastly, the ID affine coordinates associated to the nonquadrilateral (top) vertex and the perpendicularly 
projected point to the quadrilateral face are 

(9.5) 

(9.6) 

With these tools in the arsenal, it is possible to find the desired 3D aftine-like coordinates. The first 
key observation is that each vertex in the quadrilateral face is associated to four lower dimensional affine 
coordinates. The associated coordinates are those which take the value 1 at the given vertex. For example, 
vertex vi = (0,0,0) is linked to the coordinates r'o(^i)C )5 Fo(i^) and ^o(i^)- To find a 

global coordinate associated to any vertex, the idea is to combine these components such that they vanish 
at all disjoint edges and faces. One possibility is to consider the product of all four coordinates. However, 
this gives a high order function, which is somewhat inconsistent with what one would expect. Hence, the 
global coordinate should look as “simple” as possible. Fortunately, there is such a coordinate, which in fact 
has a dual interpretation with respect to its associated coordinates. It is the product of a ID scaled affine 
coordinate and the complementing 2D affine coordinate. For vertex vq, it would either be /ro( j^)t'o(C 2 , C) 
c>i' 0- These two interpretations coincide and define the pyramid affine-related coordinates. 


//o(C) = i-C, 


Fi(C) = C- 


Their gradient is 


V/io(C) = 0 , 

-1 


V/Ui(C) = 0 


78 



For the nonquadrilateral vertex, = (0, 0,1), there is an already existing affine-related eoordinate whieh 
is merely i'2{C) = Iri summary, the pyramid affine-related eoordinates are 


Ai(e,C) = Fo(i^)z^o( 6,C) = Fo(i^)i^o(6,C) = 
A2(e,C) = w(i^)^^o(6,C) = Fo(i^H(6,C) = 
A3(e,C) = w(i^)^^i(6,C) = w(i^H(6,C) = 
A4(e,C) = Fo(i^)z^i( 6,C) = Fi(i^)^o(6,C) = 

AsCC) =1^2(0 =Fi(C) =C- 


Their gradient is 


vAi(e,c) = 


/ -( 1 - 6 - 0 ' 

/ 1-c 

-( 1 - 6-0 
i-C 

\(T3^ h 


VA2(e,C)= i3 


vA4(^,C) = 


-6 
1-c 
( 1 - 6-0 
1-C 
\ -giO 
\ ( 1-02 



vA3(e,c)= ft 


va5(C) = 



(9.7) 


(9.8) 


Apart from being produets of lower dimensional affine eoordinafes, fhe pyramid affine-relafed eoordi- 
nafes fruly do behave in many ways like 3D affine eoordinafes. Firsfly, nofiee fhaf by eonsfruefion fhe fraees 
over adjaeenf faees and edges are fhe eorresponding verfex funefions of fhose lower dimensional fopologieal 
enfifies. For example, fhe fraee of Ai over faees 125 and 145 is a 2D friangle affine eoordinafe assoeiafed 
fo fhaf verfex, while fhaf of faee 1234 is a bilinear quadrilaferal vertex funelion. Seeondly, nofe fhaf every 
(^, C) in the pyramid ean be expressed as a eonvex eombinafion of fhe verfiees wifh fhe affine eoordinafes 
being fhe weighfs, 

^5 5 

To'j = ^Aa(^,C)na, wifh ^Aa(^,C) = l and 0 < Aa(?,C) < 1, (9.9) 

and where Va are fhe eoordinafes of verfex a. The main differenee wifh fhe legitimate simplex affine eoordi¬ 
nafes radieafes in fhe faef fhe fhe pyramid affine-relafed eoordinafes are nof defined by fhe properfies above 
(see §1.9). Indeed, even fhough fhey have polynomial fraees af fhe boundary, fhey involve rafional polynomi¬ 
als in fhe interior, and fhis is an inherenfly new properly. Neverlheless, for many praelieal purposes, fhey ean 
be Ihoughf of as affine eoordinafes, and from now on will be referred fo as the pyramid affine eoordinafes. 

An imporlanl remark is fhaf all fhe resulfs assoeiafed fo fhe definilions of fhe aneillary funefions were 
proved in a very general selling fhaf eneompasses fhe pyramid affine eoordinafes and fhe fael fhey ean be 
rational. In parlieular, fhe proofs of Lemmas 1 and 4 hold. 

Note fhaf all fhe affine eoordinates illuslraled ean be eompuled for pyramids wifh a parallelogram base. 
In fael, if is very easy fo make Ihese ealeulalions for pyramids wifh an arbilrarily plaeed lop vertex and whose 


79 



rectangular base is normal to the vertical ( direction and aligned with the ^ coordinates. This assertion 
includes any of the typical master pyramids found in the literature. With the affine coordinates computed, 
it is just a matter of substituting them (and their gradient) into the expressions for the shape functions to 
be presented throughout this section, so that in fact these expressions are independent of the choice of the 
master pyramid. Hopefully, this motivates other researchers to communicate their results in terms of affine 
coordinafes as well. 

Finally, by consfrucfion, fhere are nafural relationships befween fhe fopological enfifies and fhe differenl 
fypes of affine coordinafes defined. The relafed affine coordinafes are fhose which fake fhe value 1 af fhe 
prescribed topological enfify. The top vertex, v^, is linked to t'2(C) = Fi(C) = ^5(0- The quadrilateral 
verfices are each associafed to two ID scaled affine coordinafes, two 2D friangle affine coordinafes and one 
3D pyramid affine coordinafe. Meanwhile, friangle edges are linked fo two ID scaled affine coordinafes, 
while mixed edges are associated to one ID scaled affine coordinafe and fhe verfical ID affine coordinafe 
//o(C)- Lasfly, friangle faces are linked fo one ID scaled affine coordinafe, while fhe quadrilaferal face is 
linked fo fhe verfical ID affine coordinafe //o(C)- As usual, fhese associafed affine coordinafes can acf as 
nafural blending funclions. 

Exact Sequence 

If should be clear by now fhaf fhe pyramid has a fundamenlally differenl slrucfure lhan fhe previous elemenls, 
and one would expecf fhis fo have an impacf on fhe discrefe spaces fhaf affempf fo approximate fhe energy 
spaces in (1.7). 

Firsfly, nofe fhaf an absolute requiremenf is fhaf fhe frace of fhe spaces over fhe faces span fhe lower 
dimensional discrefe polynomial spaces for fhe friangle and quadrilaferal respecfively. This is whaf en¬ 
sures fhaf fhe shape funclions are compalible over adjacenl elemenls. However, any affempf al finding a 
Ihree dimensional polynomial space salisying fhose properties is futile, since one can find counterexamples 
malhemalically showing fhaf fhis lask is impossible. 

Hence, fhe use of rational polynomial spaces is fhe nexl nafural slep. This issue already arised, al 
leasl inluilively, while analyzing fhe desired properties of affine coordinafes, because fhe use of scaled 
coordinates was required. Neverlheless, dealing wilh rational polynomial spaces is difficull, and finding 
finile dimensional higher order spaces satisfying all fhe desired frace, exacl sequence and approximabilily 
properlies is a far from Irivial lask. In facl, only unlil recenlly did such conslruclions slarled fo appear in fhe 
lileralure. In fhe conlexl of fhis work, perhaps fhe besl suited sel of such spaces is fhaf proposed by Nigam 
and Phillips (2012), which is consisfenl wilh fhe “nafural” firsl order spaces analyzed firsl by Gradinaru and 
Hipfmair (1999). 

Respecling fhe nofalion of Nigam and Phillips (2012), fhe discrefe rational polynomial spaces approxi¬ 
mating (1.7) are, 

^{0),P ^(2),p ^(3),p ^ 
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where the m in eorresponds to the order of the differential form in 3D, so that the elements in 

are 0-forms, and so on. The preeise definitions of these spaees are somewhat teehnical and will be postponed 
to Appendix B. In fact, the proofs that the shape funtions lie in the desired space are also technical and 
inconveniently load the readibility of the document, so they are presented in Appendix B as well. This by no 
means implies that the spaces are not important and do not play a role in the construction. In fact, quite the 
opposite. The spaces are so well suited to the pyramid, that most of the time they impose little restrictions 
on the intuitive constructions presented here. Hence, in many ways, despite looking complicated, they are 
“natural”. 

Finally, it is worth emphasizing that the goal in this section (and in general in this work) is to motivate 
the construction of the shape functions through geometrical arguments (via the affine coordinafes defined 
before) combined wifh fhe carefully chosen ancillary operators defined fhroughouf fhe documenf. This 
approach leads fo shape functions salisfying fhe desired frace properfies and which eifher are in fhe desired 
space or can be nafurally fweaked fo lie in fhe space. The nofable excepfion is fhaf of fhe H (div) friangle 

faces, in which fhe space fruly plays a nonfrivial role and forces fo consider a more infricafe yef consisfenf 

consfrucfion. 

9.1 Shape Functions 

The dimension of fhe space is + 3p + 1. The number of linearly independenf shape funcfions will 
coincide wifh fhaf dimension. 

9.1.1 Vertices 

The verfex shape funcfions will be precisely fhe associafed 3D pyramid affine coordinafes. Indeed, lake for 
example vertex vi, so fhaf fhe vertex function is 

The frace properfies are salisfied by consfrucfion and are shown explicifly nexl, 

C)l6=o = = z^o(6,C), 

0"(e,C)la=i-c = = 0’ 

0"(e,C)l6=i-C = Ai(C,C)l^„(^)=o = = 0’ 

01^=0 = ^ , 

<f'{^,C)k =0 = Ai(.f,C)lc=o =tio(6)F(6)- 

The funclion is also in fhe lowesf order space Similar argumenfs apply to all ofher quadrilateral 

verfices and fhe lop verfex as well. 
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More generally, the vertex funetions and their gradient are, 


C) = Aa(?, C), C) = VAa(e, C), 


for a = 1, 2, 3,4, 5. There are a total of 5 vertex funetions (one for eaeh vertex). 


(9.11) 


9,1.2 Edges 


Mixed Edges. Take for example mixed edge 12. The first naive approaeh is to use the 3D pyramid 
affine eoordinafes direefly on whieh gives 


for i = 2,... ,p. This affempf almosf works beeause if is in fhe eorreef spaee, safisfies fhe vanishing 
eondifions, and even has fhe righf form af fhe edge ifself. Indeed, af friangle faees 235 and 145, i^o{P,i,C) = 0 
and 1^1 (Cl, C) = 0 respeefively, while /io(i%^) = 0 af faee 435. However, fhe nonzero fraee over fhe 
adjaeenf quadrilaferal faee is nof of fhe eorreef form, sinee if blends nonlinearly wifh fhe faelor /ro(C 2 )* 
instead of linearly like /Uo(C 2 )- Therefore fhe funelion violates dimensional hierarchy and does nof work 
for our purposes. Neverfheless fhis analysis ellucidafes how fo fix fhe issue. The idea is fo have fhe facfor 
separated as a blending factor, so fhaf fhe effecls of ^o(i^) are essentially separated from fhose 
of fi)i(Ci) C) iri Ai 2 (C, 0- Hence, fhe shape functions for fhis edge are 


C) = (i7oi(Ci, 0) = /^o(i^)(i^o(Ci, C) + to(Ci, C))* <J>i[ 


.12 


blend 


project 


evaluate 


for i = 2,... ,p. As wifh fhe previous candidafe all vanishing properties are salisfied, buf fhis lime fhe 
nonzero fraee properties are also easily seen fo hold. The projection being implied is 


( 6 ) 6)0 


( 6 , 0,0 

fAte Aa. o') 
vi-C’ 


( 


gi 

i-C 


, 0 , 0 ). 


If is a Iwo step projection, where fhe firsl step is to projecl fo an adjaeenf face and fhe second is to projecl 
along fhaf face fo fhe given edge via fhe sfandard 2D edge projecfions (see Figures 4.2 and 5.2). If fhe face 
projecfion is chosen as fhe friangle, fhen fhe projeclion af play is called fhe horizontal friangle face projection 
and consisfs of finding fhe inlerseclion P’ = (Ci, 0, 0 of the face with the projecting line parallel to the C 2 
direction and passing through the original point P = (Ci, 6, 0- This is shown in Figure 9.2. If the face 
projection is chosen as the quadrilateral, then the projection is simply the intersection P' = (|^, 0) of 

the face with the projecting line passing through the top vertex and the original point P = (Ci, C 2 , 0- This 
is shown in Figure 9.4. 
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Horizontal Triangle Face Projection 

Figure 9.2; Horizontal triangle face projection from P to P' followed by a mixed edge projection from P' to P". 


In general, the shape functions and their gradient are 

</>!(?, C) = /^c(i%)(/'f(^0l(?a,C)), 

V . (9-12) 

where i = 2,... ,p, {a, b) = (1, 2), (2,1) and c = 0,1. There are p — 1 edge function for each edge, for a 
total of 4(p — 1) mixed edge functions. 


Triangle Edges. For instance, take triangle edge 15. Again, the naive approach is to use the 3D 
pyramid affine coordinates on cpf, leading to the shape functions, 

C) = -/>f (Ai5(?, 0) = (Al(?, C) + A5(C))>f ( C) ) , 

'- 

'-V-' 

evaluate 


for j = 2,... ,p. In this case it works perfectly well, with the trace properties being satisfied. Indeed, 
Ai(?) C) = 0 over faces 235 and 435, while A 5 (^, C) = 0 over fhe quadrilateral face. Moreover fhe resfricfion 
of Ai 5 (^, C) over fhe faces 125 and 145 gives Fo2(^i,C) and Fb2(^2,C) respectively, so fhe nonzero fraces 
are fhe appropriafe friangle fraces. The projecfion being implied here is highly nonfrivial. If is a fwo sfep 
projecfion given by 




U1-S2)(1-C) 


0 , 


c 

l-?2 


0 0 __ 


The firsl step is called an oblique friangle face projecfion and consisfs of running a plane fhrough fhe original 
poinf P = (^1 , ^2 , C) and fhe opposife boffom edge fo fhe face (edge 43), followed by finding fhe intersection 
of fhis plane wifh fhe planes passing fhrough fhe ofher fwo adjacenf friangular faces (faces 235 and 415 wifh 
equations = 1 — ^ and = 0 respectively). Call fhis intersection C = (0, —1)- Finally, fhe 
infersecfion of face 125 wifh fhe projecfing line from fhe original poinf P fo fhe infersecfion C is found and 
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Oblique Triangle Face Projection 


Figure 9.3: Oblique triangle face projection from P to P' followed by a triangle edge projection from P' to P”. 

labeled as P' = 0, p:^)- This projection is illustrated in Figure 9.3. The final step is simply 

to project as usual along the 2D triangle face to the point P” = (0, 0, )■ 

In general, the shape functions and their gradient are 

mx) = ^fXaX^x)) , vmx) = V 4 >fxa 5 xx)) , ( 9 . 13 ) 

for z = 2 ,..., p and a = 1, 2,3,4. There are p — 1 edge functions for each edge, giving a total of 4(p — 1) 
triangle edge functions. 



9,1.3 Faces 




Quadrilateral Face Projection 
Figure 9.4: Quadrilateral face projection from P to P' followed by a mixed edge projection from P' to P" 


Quadrilateral Face. It was already mentioned that the quadrilateral face projection, illustrated in 
Figure 9.4, takes an arbitrary point P = {C 1 X 2 X) to the point P' = along the face. This 

projected point is actually represented by the affine coordinate quadruple (poi(p^), /roi(pir^))- Hence, the 
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natural choice is to use the quadruple with the ancillary function (/>? . This already satisfies all the necessary 
trace properties, except at the top vertex itself, where there might be a singularity. This is corrected by 
adding a factor of which also ensures the function is in the correct space. 

The shape functions and their gradient are 

4i(C,C) = , 

^ ^ (9.14) 

where f = 2,..., p and j = 2,..., p. Naturally, there are (p — 1)^ shape functions for the quadrilateral face. 


Triangle Faces. Similar to the mixed edges, there are two possibilities. Obviously they both involve 
Take for example triangle face 125. The first alternative is to use the 3D pyramid affine coordinates 
directly, yielding as a result 


C) = (Ai25(C, 0) = (Ai(C, C) + HC, C) + A5(C))*+'>; 




Ai («,C)+A2kc)+A5(C) ^125 (C, C) 


blend 


project 


evaluate 


where i > 2 and j > 1. In this case the projection implied is precisely the oblique triangle face projection 
illustrated in Figure 9.3. This function lies in the correct space and is easily seen to satisfy the necessary 
trace properties (see (5.11)). Hence, it is a perfectly valid candidate. 

A second candidate relies in the same approach taken for the mixed edges, in which the effects of the 
components of Ai 2 (C, C) separated. In that case, the functions are 

= Fo(i#c)^§( foi2(6,C) ) , 

blend project 

-V-^ 

evaluate 

for i > 2 and j > 1. Here, the projection implied is the horizontal triangle face projection shown in Figure 
9.2. Again, the function is in the correct space and satisfies fhe required trace properties, so it is also a valid 
candidate. 

The second alternative is chosen, so the general shape functions and their gradient are 

C) = Fc( (t7oi2(Ca, C)) , 

V4 (e,C) = /Uc(l^)V,/.A(Foi2(ea,C)) + 0^(^?O12(^a,C))Vpe(l^) , 

where f > 2, y > 1, n = i + j = 3,... ,p, (a, 6) = (1, 2), (2,1) and c = 0,1. As with all triangle edges, 
there are ^(p — l)(p — 2) face functions for each face, for a total of 2(p — l)(p — 2) triangle face functions. 
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9,1,4 Interior Bubbles 


The interior bubble funetions resemble elosely the ease of the hexahedron bubbles, and they are deduced 
from the quadrilateral face functions, where the factor (/roi(C)) is used instead of /Uo(C)- This ensures all 
the vanishing properties are satisfied. 

The bubble functions and their gradient are 

<t>%k (?> C) = (mi (), Ml ((mi (C)), 

V</.!(,fc(e,C)=</>f(Mi(C))V</>°(Mi(i^),Mi(i^))+</>°(Mi(i^),Mi(i^ 

where i = 2,... ,p, j = 2,... ,p and k = 2,... ,p. Clearly there is a total of {p — 1)^ interior bubble 
functions. 


(9.16) 


9,2 i?(curl) Shape Functions 

The dimension of the space is 3p^ + 5p. The same number of shape functions will span the space. 

The construction of the shape functions for H (curl) is completely parallel to that of and they involve 

the same underlying projections. 


9,2,1 H (curl) Edges 

iJ(curl) Mixed Edges, Take for example mixed edge 12. As in using m(i^) und f(Di(Ci)C) 
separately instead of Ai 2 (^, C)? it follows the edge shape functions are 

^1(^,0 = M(i#c)^f(^oi(a,c)) = M(i%)(^o(6,c)+z^i(6,c)r+2i^f 


blend 


1^0 ( 6 , 0 + 1^1 C ) 

'-V-' 

project 

V 

evaluate 


for f = 0,..., p — 1. From the H (curl) edge triangle functions, it follows that along the edges 15 and 25 the 
tangential component of Ef{uoi{S^i, Q) vanishes, and due to its independence from ,^2 it immediately fol¬ 
lows that the same is true for the faces 235 and 145. Along face 435, it holds that M(i^) = 0 so that it also 


vanishes there. Along face 125 /ro(i^) = und the function becomes (Fbi(^i, C)), which as desired 


IS 


C 

I, 

the triangle 2D trace for the edge functions. Finally, at the quadrilateral face, po(i%^) = ^- 0 (^ 2 ), while the 
tangential component of £'^^(^ 01 (^ 1 , C)) is the corresponding segment ID edge function, meaning that 
the trace along this face is iJ.o{C2)Ef{j2oi{^i)), as required. Hence, all trace properties hold. The projection 
implied is the mixed edge projection depicted in Figures 9.2 and 9.4. Lastly, the shape functions are in the 
correct space. 

The shape functions and their curl are 

V X £|(C,C) = m(i^)V X £f (Foi(ea,C)) + Vp,(^) X , 


(9.17) 
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where i = 0,... ,p — 1, {a,b) = (1, 2), (2,1) and c = 0,1. There are p edge funetions for eaeh edge, for a 
total of Ap mixed edge functions. 


iT(curl) Triangle Edges. For instance, take triangle edge 15. Like in H^, one can directly use the 3D 
pyramid affine coordinates on Ef. The resulting shape functions are 


E^(C,C) = (Ai5(e,C)) = (Ai(e,C) + A5(C))*+'K’' 


blend 


Ai(«,C)+A5«)^15(^,C) 

'-V-' 

project 


evaluate 


for i = 0,..., p — 1. To argue the nonzero traces have the correct form, take for example face 125, and the 
decoupling Ai(.^, C) = 0- Then, using that t'2(C) = A 5 (C). the lowest order element is 

^0^(Ai 5(C,C)) = /Uo(i#j)i^o(ei,C)VA5(C) - A5(C)v(po(i^)z^o(ei,C)) 

= (^02(6,0) - ^^0(6, C)z^2(C)Vpo(i^) • 

When evaluated at face 125, as a result Vpo(i^) is orthogonal to the face (it is an isosur¬ 

face), so that the tangential component at the face is precisely the nonzero components of i?^(Fo2(^i, C))- 
Hence, its nonzero trace on the face is a triangle 2D H (curl) edge function, as expected. Using the same 
argument but at face 435, where = 0, this time the tangential component vanishes completely. 

Symmetric arguments apply to faces 145 and 235. Finally, at the quadrilateral face, where A 5 (C) = 0, 
the lowest order function takes the form Ai(^, C)VA 5 (C) which is normal to the face (it is an isosurface of 
AsCC)), so the tangential component is zero. These arguments confirm fhaf fhe frace properties hold. Lasfly, 
fhe projecfion is a friangle edge projecfion as illusfraled in Figure 9.3. 

The shape functions and fheir curl are 


EliCX) = EfMtO) 


VxE!itC) = VxEfiXa,{(,C)), 


(9.18) 


for i = 0,... ,p — 1 and a = 1, 2, 3,4. There are p edge functions for each edge, for a fofal of 4p friangle 
edge funcfions. 


9.2.2 H (curl) Faces 

H (curl) Quadrilateral Face. As expecfed, fhe projecfion implied in fhese expressions will be fhe same 
as fhaf of If is a quadrilaferal face projecfion as depicfed in Figure 9.4. However, in fhis case fhere 
will be fwo families. Bofh families will easily satisfy fhe vanishing frace properties by use of (3.7), (4.10) 
and fhaf Vpc(j%^) is orfhogonal fo fhe friangle faces where = 1- The only difference wifh 

radicates in fhe use of fhe higher order blending function pq{C)‘^, which is used in order fo be in fhe correcf 
space. There is a grand fofal of 2p{p — 1) quadrilaferal face functions. 
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Family I: The shape funetions and their eurl are 


4(^,0 = , 

V X 4(^,0 = Mo(C)'V X (9.19) 

+ 2Mo(C)V/ro(C) X , 

for i = 0,...,p — 1 and j = 2,... ,p. There are p(p — 1) shape functions in this family. 

Family II: The shape functions and their curl are 

V X Elj{^,C) = po(C)'V X (9.20) 

+ 2po(C)Vpo(C) X i?°(poi(i%^),m(i^)), 

for i = 0,...,p — 1 and j = 2,... ,p. Note the fact that the entries are permuted with respect to the first 
family. There are p(p — 1) shape functions in this family. 

/^(curl) Triangle Faces. As with there will be two valid alternatives. Take for example face 125. 
The first alternative is to use the pyramid affine coordinates directly on Efj and multiply by po(i 3 ^)~^. This 
approach is discarded in favor of separating the effects of po(i^) and (^) directly from Ai 2 (^, C). 

The resulting functions satisfy the necessary vanishing conditions using similar arguments to those used 
for mixed and triangle edges. Also, the projection implied is the horizontal triangle face projection shown 
in Figure 9.2. There are p(p — 1) functions per triangle face, for a grand total of 4p(p — 1) triangle face 
functions. 


Family I: The shape functions and their curl are 


V X 4(^,0 = Pc(i^)V X Ef^ii;ou{^a,C))+yd^ci^c 

for f > 0, j > 1, n = i + j = 1,... ,p — 1, (a, 6) = (1, 2), (2,1), and c 
functions in this family. 


A 

) X Efj{i7oi2{(a,C)) , 

= 0,1. Every face has ^p(p — 1) 


Family II: The shape functions and their curl are 

0 = Pc{^^)E^{dMCa, 0) , 

V X 4(^,0 = Pc(i^)V X + Vpe(i#^) X Ef^XM^aX)), 

for i > 0, J > 1, n = f + j = 1,.. . ,p — 1, (a, 6) = (1, 2), (2,1), and c = 0,1. Note the fact that the entries 
are duoiCa, C) as opposed to Foi 2 (^a, C). Every face has ^p(p — 1) functions in this family. 



9.2.3 H (curl) Interior Bubbles 

These will be separated aeeording to a Helmholtz deeomposition. Indeed, there are four families of interior 
bubbles, and the first family eorresponds preeisely to the gradients of interior bubble funetions, so they 
have zero curl. The other three families will essentially be generated by quadrilateral face functions in 
H (curl) and In all cases, the trace properties follow easily. There is a grand total of 3p{p — 1)^ interior 

bubble functions. 

Family I: These are the gradients of interior functions. The shape functions and their curl are 

= </>f(m(C))V0°(m(i^),m(^)) 

+ </>° (mi (), /2oi ((/2oi (C)), (9-23) 

Vxi^b,(e,C) = 0, 

where i = 2,... ,p, j = 2,... ,p and k = 2,... ,p. There are {p — 1)^ interior bubble functions in this 
family. 

Family II: The shape functions and their curl are 

=M(C)0f(Mi(C))^°(Mi(i^),Mi(^^ , 

V X = M(C)0f(Mi(C))V X i?°(^oi(i^),Mi(i%c)) (9-24) 

+ (Mo(C)V</)f(/2oi(C)) + </*f(Mi(C))VM(C)) X , 

for i = 0,...,p — 1, j = 2,... ,p, and k = 2,... ,p. There are p{p — 1)^ shape functions in this family. 
Family III: The shape functions and their curl are 

EfjkXX) = MC)4>kimiC))E?j(^fioii^),miSc)) ’ 

V X = M(C)0f(Mi(C))V X i?°(^oi(i%^),Mi(i^)) (9.25) 

+ (po(C)V</)f(/2oi(C)) + </>f(Mi(C))VM(C)) X i?°(Mi(i#f),Mi(i^)) , 

for i = 0,... ,p — 1, j = 2,... ,p, and k = 2,... ,p. Note the entries are permuted with respect to the 
second family of interior bubbles. There are p{p — 1)^ shape functions in this family. 


Family IV: The shape functions for the final family of interior bubbles and their curl are 

(Mi(i%^),Mi(i^))nM(C)"-'V/io(C), 

V X F;’)(C,C) = nM(C)"“'V,/.°(poi(i^),Mi(i^)) x Vpo(C), 
for i = 2,... ,p, j = 2,... ,p, and n = max{i, j}. There are {p — 1)^ shape functions in this family. 


(9.26) 
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9.3 i?(div) Shape Functions 

The dimension of the spaee ^ is 3p^ + 2p. The number of linearly independent shape funetions eoineides 

with that dimension. 

9,3.1 H (div) Faces 

iT(div) Quadrilateral Face. These are eonstrueted exaetly the same way as the and Tf(curl) eoun- 
terparts, but with the higher order blending funetion po{C)^ so that the resulting funetions lie in the spaee. 
In view of (6.14), it is elear that the funetion will only have a tangential eomponent along the triangle faees, 
so that the normal traee vanishes at these faees, as required. Moreover, again through (6.14), it is evident 
that the nonzero traee on the faee is a quadrilateral faee funetion. The projeetion is onee again depieted 
in Figure 9.4. 

In view of (6.14), the shape funetions and their divergenee are 

F^(e,C) = Fo(C)V°(;uoi(i^),m(i^)), 

, ^ , (9.27) 

V • F^(e,C) = 3/ro(C)'VMo(C) • , 

where i = 0,...,p — 1 and j = 0,... ,p — 1. Clearly, there are quadrilateral faee funetions. 

H{dW) Triangle Faces. The eonstruetion of the T7(div) triangle faee funetions is highly nontrivial. 
Indeed, an analogous eonstruetion to the ease of or H (curl) does not work here. This radieates in the 
definition of the spaee itself. The issue is even present for the lowest order spaee, and in faet it is by looking 
at this spaee in detail that the problem is solved. 

To summarize the eonstruetion, take for example faee 125. The detailed ealeulations are in Appendix B. 
Proeeeding as in the or H (curl) case, one obtains the following two disheartening facts for the lowest 
order candidate functions, 

Fo(i#c)^oo(^012(6,0) Fo(i^)“Vo^(Ai25(e,C)) 

However, they both satisfy the desired trace properties. One could alternatively attempt a more direct con¬ 
struction, but issues arise constantly, either because the functions are not in the space, or because there are 
“illegal” derivatives of Legendre polynomials Pi, which in theory should not exist as they are elements of 
and are intended to approximate elements of that space (for example, a discontinuous function). These 
issues do not arise when using in view of Lemma 4, and this is one of the reasons why it is so conve¬ 
nient to use it. Fortunately, there is a way to make this happen. The key is to look at the unique lowest order 
function for a given face. The explicit formulas are given in Gradinaru and Hiptmair (1999) and Nigam and 
Phillips (2012). After scrupulous observation, one obtains 

^oo(C,C) = ^(Fo(i^)Fo^(i7oi2(6,C)) + Fo(i^)“Vo^(Ai25(C,C))) 
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Hence, this suggests the following general formula for the face 125 shape functions, 

VU^X) = ^(f^0(i%^)V^.f(l7012(6,C)) + M0(l^)“VA(Ai25(^,C))) , 

for i > 0, j > 0, and n = i + j = 0,... ,p — 1. In Appendix B it is shown that these high order 
functions are in the correct space and that they sastisfy the trace properties. Some could worry when seeing 
the factor but this is in fact not a real singularity. Indeed it is shown in Appendix B how to 

avoid it explicitly, along with an alternative formula convenient for computations. Lastly, note the inherent 
projection is not unique, and in fact is a combination of horizontal and oblique triangle face projections. 

Finally, in view of (7.14) the general shape functions are 

V-F,f(C,C) = •^,f(%2(ea,C)) + Fc(i^)“V-H,f(A,e5(C,C)) (9-28) 

- • FA(A,e5(C,C))) , 

where i > 0, y > 0, n = i + j = 0,... ,p — 1, (a, 6) = (1, 2), (2,1), c = 0,1 and where (d, e) depends on 
(a, b, c) (in fact \debX-, C) = Fc(i^) j^oi2(Ca, C))- There are \p{p + 1) functions for each face, leading to 
a total of 2p(p — 1) triangle face functions. 

9.3.2 H (div) Interior Bubbles 

Like the H (curl) interior bubbles, these will be separated according to a Helmholtz decomposition. Indeed, 
the first three out of seven families of interior bubbles are precisely the curl of H (curl) interior bubble 
functions, so they have zero divergence. In all cases, the trace properties follow easily. There is a grand total 
of 3p‘^{p — 1) interior bubble functions. 

Family I: The shape functions and their divergence are 

vU^X) = Fo(C)0f(m(C))v X 

+ (Fo(C)v</.f(m(C)) + 0f(m(C))VMo(C)) (9-29) 

v-v^^,XX) = o, 

for f = 0,...,p — 1, j = 2,...,p, and k = 2,... ,p. There are p(p — 1)^ shape functions in this family. 

Family II: The shape functions and their divergence are 

v^.kXX) = Fo(C)</>f(m(C))v X 

+ (Fo(C)v</.f(m(C)) + 0f(m(C))Vpo(C)) (9-30) 

v-Fb,(e,C) = o, 
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for i = 0,... — 1, j = 2,... ,p, and k = 2,... ,p. Note the entries are permuted with respeet to the first 

family of interior bubbles. There are p(p — 1)^ shape funetions in this family. 


Family III: The shape funetions and their divergenee are 

= n/io(Cr“'V</.° (;U0i(i^),m(i^)) X Vmo(C) , 

V-F,^(C,C) = 0, 

for i = 2,... ,p, j = 2,... ,p, and n = max{f, j}. There are {p — 1)^ shape funetions in this family. 


(9.31) 


Family IV: These have nonzero divergenee and are generated by the quadrilateral faee funetions, but 
using ;Uo(C)^0f (Moi(C)) a faetor instead of pq{C)^. In view of (6.14), the shape funetions and their 
divergenee are 

v-F,^fc(^,C) = (/uo(C)'v</.f(^oi(C))+2fro(C)</>f(m(C))VMo(C))-^°(m(i^),m(i^)), 
for i = 0,...,p — 1, y = 0,...,p — 1, and k = 2,... ,p. There arep^(p — 1) shape funetions in this family. 


(9.32) 


(9.33) 


Family V: These have nonzero divergenee and are expressed as the produet of a power of pi{C) with 
a eurl. As a first step, define 

= (So^i)) +foVfox((/>f (so"i)V(/>f , 

for i = 2,... ,p, and j = 2,... ,p. Clearly, V • V-j{sqi, Sg^, fo) = 0. The shape funetions and their 
divergenee are 

V'.'e.c) = m(C)”-'v'j(wi(A:).*i(A).'‘o(C)), 

(9.34) 

V • v;5(c,c) = (n - i)pi(c)"-'Vpi(C) • (poi(i^),m(i#j),//o(C)), 

for f = 2,... ,p, y = 2,... ,p, and with n = max{i, j}. There are (p — 1)^ shape funetions in this family. 


Family VI: These have nonzero divergenee and are expressed as the produet of a power of pi(C) with 
a eurl. First define 

^i“(soi,Fi,io) = v(^fo(/>f(soi)) X Vpi = (^foV(/>f(soi) + 2fo</>f(soi)Vfo) x Vpi, (9.35) 

for f = 2,... ,p. Obviously, V • V^“(soi) Ft) ^o) = 0. The shape funelions and their divergenee are 

V ^ ^>(5,0 = (i - i)A.i(cr"v,.i(C) ■ rf (wi{^),,.i(A;). w(C)). 

for f = 2,..., p. There is a total of (p — 1) shape funetions in this family. 


(9.36) 
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Family VII: Using (9.35), the shape functions and their divergence are 

^ ^ (9 37) 

V • C) = (j - i)/ri(cy-'Vw(C) • m(i^),f^o(C)), 

for j = 2,... ,p. Note the entries are permuted with respect to the sixth family of interior bubbles. There is 
a total of (p — 1) shape functions in this family. 

9.4 Shape Functions 

The dimension of the space is p^. The same number of shape functions will span the space. 

9.4.1 Interior 

Again, these are reminiscent of the shape functions for the hexahedron. They are, 

V’?,fc(e,C) = [^d(m(i^))[^’i](m(i^))[^’A:](m(C))(Vi^i(6,C)xVi.i(6,C))-Vw (9.38) 

for i = 0,..., p — 1, j = 0,..., p — 1 and A: = 0,..., p — 1. There are p^ interior functions. 

9.5 Orientations 



Figure 9.5; Master pyramid with numbered vertices and local edge and face orientations. 

The predefined local edge and face orientations for the pyramid are illustrated in Figure 9.5. They 
represent the o = 0 case. The task at hand is to find the associated locally ordered tuples of affine coordinates 
representing those local orientations. The key is being aware of the relationships between the vertices and 
the affine coordinates, which have been explained in detail at the beginning of this section. Hence, once §6.5, 
§7.5 and §8.5 are consulted, it should be clear how to construct the orientation embedded shape functions. 
Consult Appendix B if the reader wants to avoid computational instabilities for the H (div) triangle face 
functions. 
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10 Conclusions 


We have presented here a full systematie eonstruetion of hierarehieal higher order shape funetions for ele¬ 
ments of “all shapes” (see Figure 1.1) using the exaet sequenee logie. Compatibility of the shape funetions 
at the interelement boundaries is based on the idea of having a known fixed traee at the boundary whieh 
is extended (or lifted) to the rest of the element. Henee, the shape funetions ean be used in hybrid meshes 
eontaining elements of all shapes. Furthermore, due to the properties of the diserete spaees, interpolation 
estimates are ensured in any hybrid master element mesh and for all energy spaees. 

The unified eonsfruefion is based af ifs eore in eonsidering fensor produefs of polynomials. This has pos- 
ifive implieafions from fhe eompufafional poinf of view, and eould resulf in fhe sueeessful implemenfafion 
of fasl integration feehniques as deseribed in Appendix C. 

Also, fhe shape funelions allow fhe polynomial order p fo vary aeeross a given mesh. For example, fhe 
quadrilaferal, hexahedron and prism shape funelions are nalurally anisolropie and can have differenl orders 
in each direclion. More so, each edge and face in fhe mesh can have Iheir own order p, independenl of fhe 
order of fhe neighboring edges, faces and inferiors in fhe mesh. Hence, feehniques lhal exploif fhe use of 
local p adapfivify can be implemenfed using fhese shape funelions. 

As polynomial building blocks, Legendre and Jacobi polynomials are used in Ihis work. Our resulls 
show lhal fhe recursive formulas presenled in Ihis documenl lo implemenl Ihose polynomials seem lo be fhe 
mosl accurate in comparison lo olher recursive formulas. The choices of Jacobi and Legendre polynomials 
are known lo have exlremely good sparsily and conditioning properties for typical projection problems 
(Beuchler el ah, 2012a). However, Ihere is flexibility in Ihis choice of polynomials, and il mighl be worlh 
investigating if Ihere are applications where differenl choices provide useful advanlages (see Appendix A). 

All conslruclions are written in terms of affine coordinates and Iheir gradient Indeed, Ihe shape functions 
are valid for any (typical) master elemenl geomelry, provided Ihe affine coordinates are computed. We hope 
Ihis has convinced Ihe reader lhal Ihe direcl use of affine coordinates for all elemenls (nol jusl simplices) 
is Ihe ideal approach, and lhal il motivates olher researchers lo also communicate in Ihose general terms. 
The polynomials we chose are shifted to have the domain (0,1) instead of the typical (—1,1). We claim 
this is the natural choice for construction of shape functions, since all affine cordinates have range [0,1] 
and affine coordinates are the natural inputs for the polynomials. Hence, we encourage the implementation 
of the polynomials in the shifted domain. The concept of polynomial homogenization was introduced and 
heavily used. It is a particular form of scaling which is closely related to affine coordinates, and is a tool that 
provides natural extensions. In fact, homogenization provides some level of geometrical intuition, since the 
projected affine coordinates arise naturally through this process. Theoretically it is also a convenient tool 
since it results in homogeneous polynomials, which have many desirable properties. 

Moreover, only eight ancillary operators effectively generate all shape functions. These ancillary oper¬ 
ators are coordinate free, in the sense that the form of the operators is invariant with respect to any transfor¬ 
mation. This is important, because it allows to transform nonlinearly to other geometries. Hence, it suffices 
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to compute the affine coordinates and their gradient in that deformed space. Then, the shape functions re¬ 
sulting from the substitution of the deformed affine coordinafes will precisely be fhe well defined pullback 
of fhe original shape functions. This has bofh fheorefical and practical implicafions. For example, curved 
physical elemenfs are deformalions of fhe master elemenf domains. This has fhe pofenfial fo resulf in more 
efficienl compufafions when infegrafing (see Appendix C for fhe basics of integration). 

For fhe face and edge shape functions, fhe logic of projecfing, evaluafing and blending is used consis- 
fenfly for all elemenfs and all energy spaces. This provides a firm geomefrical infuifion of fhe expressions 
and formulas for fhe shape funclions, which we hope fhe reader will appreciafe. 

Addifionally, fhe shape functions can be converted fo orienfafion embedded shape functions via only 
three local-fo-global permufafion funclions (one for edges, one for Iriangle faces, and one for quadrilateral 
faces). These orienfafion embedded shape funclions are exlremely praclical in many applications, especially 
in fhe implemenlalion of conslrained nodes in hp melhods. 

All fhe characlerislics above prove fo be vilal in fhe implemenlalion of a code. Indeed, fhe number of 
imporlanl routines which are called repealedly is very small, and Ibis minimizes fhe sources of errors, while 
allowing a very focused oplimizalion of fhe implemenlalion. A complele Forlran 90 code supplemenls 
Ihis documenl.^ If provides an excellenl guidance if fhe reader is ever inleresled in implemenling Ibis 
conslruclion. The code has been lesled Ihoroughly by numerically checking polynomial reproducibilily and 
exacl sequence properlies. This is described in Appendix D. The shape functions for all elemenfs and all 
spaces are convenienlly summarized in Appendix E. 

Easily, special allenlion is given fo fhe successful conslruclion of fhe pyramid shape functions, which is 
rare in fhe lileralure. A Ihorough geomelric infuifion for fhe pyramid was described, and fhe 3D pyramid 
affine-relaled coordinafes were defined and analyzed. The sef of exacl sequence spaces were laken from 
Nigam and Phillips (2012), and Ihey are consislenf wilh fhe fundamenlal firsl order elemenfs described by 
Gradinaru and Hiplmair (1999). We believe if fo be fhe firsl lime lhal an arbilrary high order conslruclion 
of H (div) shape funclions has been implemenled whilsl respecling Ihose lower order spaces (Ihere have 
been olhers where eilher fhe lower order spaces have been larger, or simply differenl). This may prove fo 
be very valuable fo olher researchers even if only for comparison purposes. Eor fhe pyramid, if mighl be 
possible fo invesligale heller choices of polynomials for Ihe bubbles, as Ihis may provide belter conditioning 
and sparsily properties. 

To finalize, we hope Ibis conslruclion has been useful in ils melhoodology and lhal il motivates furlher 
research in Ibis very rich area. 

Acknowledgements. The work of Euenles, Keilh, Demkowicz and Nagaraj was supported wilh granls by 
AEOSR (EA9550-12-1-0484), NSE (DMS-1418822) and Sandia National Eaboralories (1536119). 


’See the ESEAS library available at https://github.com/libESEAS/ESEAS. 
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A Polynomial Families 


Let V^{x) = span{x-^ ; j = 0,..., i} be the polynomials of order i, and Fp = {Pi :i = 0, 1} and 

Fp = {P“ : j = 0,... ,p — 1} be sets of polynomials with domain x G [0,1]. The families Fp and Fp 
should satisfy that span(J^p) = span(J^") = Pi G V^{x) and Pj" G V^{x). Moreover, the family 

Fp should satisfy the zero average property, Pi(x) dx = 0, for alH > 1. 

These sets of conditions are very easy to satisfy. For example consider any sequence of polynomials 
of increasing order, say 1 ^ X ^ Xi ^ X j • • t j X P Then it is only a matter of adding a suitable constant to 
all polynomials of order i > 1 such that their integral is 0. For example, 1, x, x^, x^,..., x^ becomes 
1, X — 5 , , x^“^ — and this latter one is already a suitable family Fp. 

The elements of Fp and Fp are thought of being elements of L^, even though they are infinitely differen¬ 
tiable as polynomials. Indeed, it is often desirable that the elements of Fp and Fp satisfy certain properties 
in the (or weighted inner product, since they can result in considerably sparser finite element ma¬ 
trices. In fact, orthogonality of the Pi, is generally seeked. If there is no weight, orthogonality of the Pi 
is attained uniquely (up to a constant) by the Legendre polynomials, and this is why they are the typical 
choice. Meanwhile, the family Fp can also be chosen wisely by taking into account a weighted space 
relevant to the triangle element. If the Pi are Legendre polynomials, the natural choice for the is to be 
Jacobi polynomials with certain weights. 

Now, dehne the polynomials L* G 'P*(x) and L" G V^{x) from the Pj_i G Fp and Pj"_i G Fp as 

Lj(x) = J Pi-i{x)dx, L“(x) = J P“_i(x)dx, (A.l) 

for i = 2,... and j = 1,... ,p. Clearly, it follows that pP{x) = 'P^(x) U {Li : i = 2,... ,p} and 
pP{x) = V^{x) U {L" : j = 1,... ,p}. By construction, the Li and L“ are elements of and as a result 
their pointwise evaluation exists. In fact, due to the zero average property. 


Li{0) = Li{l) = 0, fori >2, 

L" (0) = 0 , for i > 1. 


(A.2) 


Now, apply the definition of scaling given in (2.6) to the polynomials, yielding Pi{x;t), Pf{x;t), 
Li{x; t) and L"(x; t), where x G [0, t]. In particular, note that 

Li+i(x;f) = y Pi{x)dx = fj Pi(J^^ dy . (A.3) 

Since the scaled polynomials are now bivariate polynomials in (x, t), it is useful to find fhe derivafives in 
bofh of fhese variables. Using fhe latter equality above, 

—Li+i(x;f) = = Pi{x]t). (A.4) 
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Meanwhile, using the other equality, it follows 

- (f)p(l)) . (A.« 

This suggests the definition 

Ri{x) = {i + l)Li+i{x) - xPi{x), (A.6) 

for i > 1. Using the leading order term of Pi, the reader may observe that Ri is an order i polynomial {not 
order i + 1), so the indexing with i makes sense. More importantly, 

^Li+i{x;t) = f(^{i + - (l)-Pi(l)) = Ri{x',t). (A.7) 

Exaetly the same analysis applies to the meaning that 

—Lj+ia{x;t) = P^{x;t), —L"+i(x;f) = , (A.8) 

for j > 0, and where 

Rj{x) = (j + l)L"+i(x) - xP“(x). (A.9) 

As observed, all these properties are dedueed at a very general level, and apply to any families Pp and 
Pp satisfying the simple set of properties mentioned initially. To eonelude, it follows all the results dedueed 
throughout the doeument still hold for these more general polynomial families, ineluding the vanishing 
properties, the aneillary funetion properties, ete. Henee, the reader may deeide to ehange these families if 
desired. This eould potentially provide better sparsity properties depending on the problem. Nevertheless, 
be aware that for the elassieal Laplaee problem (and many related sets of problems), the ideal polynomial 
families are the Legendre and Jaeobi polynomials, whieh are used throughout this doeument and are eonve- 
niently defined fhrough reeursive formulas. 
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B Pyramid Supplement 



Figure B.l; Master pyramid with numbered vertices. 

The master pyramid is shown in Figure B.l in the C) space {not the (Ci)^2)C) space). More 
specifically, the definition is O = {(^, ry, G : ^ > 0, r/ > 0, C > 0, ^ + C < 1, ?? + C < !}• 

This supplement provides the proofs that the pyramid shape functions proposed in §9 are in the finite 
element spaces defined in the fundamental work by Nigam and Phillips (2012). Recall from (9.10) that the 
discrete spaces forming the exact sequence are for m = 0,1, 2, 3, where m stands for the differential 

m-forms lying in each space. The spaces are defined as 

n , (B.l) 

where the spaces are called the underlying spaces, and the spaces are called the compatibility 

spaces. These two families of spaces will be defined next. 

First, note that the compatibility spaces ensure the elements in are compatible at the level of 

spaces with the other elements. Indeed, they consist of those functions having their face traces lying on 
the appropriate 2D quadrilateral and triangle spaces in (4.2) and (5.3) respectively. More specifically, let 
tr^ and tr^ be the trace of the differential m-forms over the four pyramid triangle faces and the pyramid 
quadrilateral face respectively. Recall that for m = 0 the trace is the value of the function itself, for m = 1 
the trace is the 2D tangential component, while for m = 2 the trace is the normal component. With this in 
mind, the compatibility spaces are,^ 

r(i).n = {Ee H{cnr\) : trW(i?) E trW(AAP(C, ry, C)), trg)(i?) E 

= {y E i/(div) : tr^2)(y) E trg)(7^rP(e,ry,C)), trg)(l/) E , 

r(3)-P = {r/; E L^} . 

*Note that in (B.2), the very special property that V^, and TZT^ are affine invariant is used. If the spaces were different, one 
would need to take the 2D pullback of each triangle trace to the master triangle. The same holds for the quadrilateral trace, where 
in this case it is exploited that the quadrilateral face is the master quadrilateral and no transformation is needed. 
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Fortunately, at the level of shape functions, dealing with the compatibility spaces is not as intimidating 
as it might look. All that is required is that the shape functions satisfy the dimensional hierarchy, so that 
their nonzero face traces correspond to the lower dimensional shape functions, which are known to lie in the 
correct space. Therefore, if a shape function satisfies the required trace properties, it automatically belongs 
to the appropriate compatibility space. 

The main difficulty lies in showing that the shape functions belong to the underlying spaces 
In fact, these spaces are not nicely defined direcfly on the master pyramid. For this reason, it is more 
convenient to define them on a deformed space where the symmetries are more evident, and then use the 
inverse pullback to the master pyramid. This process is explained in what follows. 



Figure B.2: Transformation from the infinite pyramid to the master pyramid. 


The deformed space is usually chosen as a cube, but Nigam and Phillips (2012) chose it to be an infinite 
pyramid, which is defined as floo = {(x, y, z) G : 0 < x < 1, 0 < y < 1, z > 0}. The first step is to 
consider the transformation r : floo —^ from the intinite pyramid to the pyramid 12, which is given by 
the component equations. 




X 


V = 


y 


C = 


(B.3) 


1 + z' ' 1 + z' ’ 1 + z 

Clearly, r is a diffeomorphism between these open sets {not when including the boundary). The inverse 


^-1 


12 —>• 12oo is given by the component equations. 


X = 


1-C 


V = 


1-C’ 


z = 


c 


1-C 


(B.4) 


The transformation is depicted in Figure B.2. Take note of the following useful expressions resulting from 
this transformation. 


1 — X = 


1-C-C 
1-C ’ 


l-y = 


i-??-C 

1-C 


1 + z 


1 

W’ 


1 

1 + z 


= !-(■ 


Next consider the following isomorphic mappings, 

^ y(m),p _ ^y(m),p'j ^ 

.y{m),p _, y{m),p 

I . y T y QQ , 


(B.5) 
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where r and are the inverse pullbaek and pullbaek mappings induced by r. Since the spaces 

and are isomorphic, it is mathematically irrelevant which of the two spaces is actually defined, 

since the other space can be determined through the corresponfing pullback mapping. However, sometimes 
there are practical reasons to explicitly define one sef of spaces over fhe ofher. Indeed, if is more opporfune 
fo define fhe deformed underlying spaces V(^ , which are 

V™'” = {•!>€ Qj;-'"-' X X ') , 

^ vxE6QJE‘’’-‘xq';5-'''»-'xc^;:’''-‘''-’ 

where = Qn'^’^ix, y, z) are fhe n-weighfed tensor polynomial spaces. They are defined as 




(B.6) 


QP’' 


’'?’’’(x,y,2:) = { G y, z)| . 


(B.7) 


Nofe fhe useful inclusion Qn'^’' C which holds for fhese rafional polynomial spaces. 

Lasfly, nofe fhaf fhe pullbacks lake differenl forms depending on m. Indeed, if Jr is fhe Jacobian 

of fhe Iransformafion r, fhe pullbacks are. 


(B.8) 


T*’^^'^(p = (po T , 

(P£H\ 

= jJ{Eot), 

E £ Tf(curl) 

T<^)v = det{Jr)JrHy ° r), 

V £ ff(div) , 

T*'^^^'ip = det(JT-)('0 o r), 

Ip £ . 


The same relafions hold for fhe inverse pullbacks buf replacing r by r“^, Jr by and fhe do¬ 
mains by fheir isomefrically isomorphic counterparls,^° iJ(curl)oo = (curl)), 

if(div)oo = {div)), and Even fhough fhey will be unnecessary, in fhe inferesf of 

complefeness fhese pullback mappings are wriffen explicilly below. 




0 




T*A‘i)V 


(pOT, 


{1 + zf 


1 


(1 + Z)3 

1 

(TT^ 


/l + z 0 
0 l + z 
\ -X -y 

/l 0 X \ 

0 1 2 / 

\0 0 l + z) 

{'ipor), 


0 I (£^or), 


(Hot), 


E £ H (curl), 

V € ff(div}, 

-tp £ . 


(B.9) 


®Note there is misprint in Nigam and Phillips (2012) in (3.8) when presenting the equivalent characterization of The 

definition here corrects that, and is consistent with the calculations in §3.3 of Nigam and Phillips (2012). 

'®The spaces are isometrically isomorphic if the appropriate weights are added to the definition of the norm. 
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The inverse pullbaeks explieitly are, 


T 


*:(i) 




(1-0^ 


1 


1 

(1^ 


/i - c 0 o\ 

0 1 - c 0 (i? O T ~^), 

V ^ ^ 1 / 

A 0 -e \ 

0 1 -7? (y o T~^) , 

\o 0 1-c/ 

(V’or"^), 


E G F(curl)oo , 

F G i7(div)oo , 


(B.IO) 


To prove the shape funetions defined in §9 are in the underlying spaees, in general one would pull 
them baek (using (B.9)) and eheek if they belong to the spaees . However, due to the coordinate 

free definitions of the aneillary funetions and the shape funetions in general, it is not neeessary to find the 
pullbaek explieitly. Instead, one simply finds the trivial pullbaek of all the sets of affine eoordinates and 
their gradients in the (x, y, z) eoordinates. Then, the only task is to evaluate all the shape funetions with 
these affine eoordinates and the result will be the desired pullbaek of the original shape funetion. This is 
why the mappings (B.9) and (B.IO) beeome redundant for our shape funetions. 

In view of these eomments, it is useful to have all the affine eoordinates and their gradients in the (x, y, z) 
spaee. The triangle affine eoordinates (see (9.1)) are 


^o{x,z) = ^ , 
^ , 

Their gradient (see (9.2)) is 


, 

i'i(y,z) = , 


1'2{Z) = , 

Mz) = . 


/-(l+^)\ /(l+z)\ 

Vl'o(x,z) =^0 ^J, Viyi(x,z) = o J, Vo 2 iz) 

Vr'o(2/,^) = , Vz/i(y, z) = ^ (i+j:) ^ , Vu2iz) 


The sets of quadrilateral sealed ID affine eoordinates (see (9.3)) are 


(B.ll) 



(B.12) 


yoix) = 1 - X, 

My) = 1 - y, 


Their gradient (see (9.4)) is 

vMx) 

'^My) 



Mx) = x, 

My) = y- 


V^i(x) 

VMi(y) 



(B.13) 


(B.14) 
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The vertical ID affine coordinates (see (9.5)) are 


Their gradient (see (9.6)) is 

= (1+^)2 = (1+2)2 ^0^ • 

Finally, the pyramid 3D affine coordinates (see (9.7)) are 

\i{x,y,z) = , A 2 (x,y,z) = , M{x,y,z) = ^, 

A4(x,y,2;) = , ^b{x,y,z) = 

Their gradient (see (9.8)) is 


V\i{x,y,z) = 


/ -(t-v) 

/ l+z 

l-\-Z 

- (1+^)2 


, VX2{x,y,z) = 


I 1 + 2 

I , '^h{x,y,z) = 


(T+ip 



/ ^ 

/ 1 + 2 


VA4(x,2/,2;) = 


^ I > VA5(z) = 
(i+P" 


/ 0 

0 

(T^ 


(B.15) 


(B.16) 


(B.17) 


(B.18) 


Before beginning with the proofs, take note of the following useful results. 

Lemma 5. Let Li and L“ be the integrated Legendre and Jacobi polynomials of order i > 2 and j > 1. 
Then there exist homogeneous polynomials [V’i- 2 ] £ si) and [Xj-i] ^ ti) such thaf^ 

[Li](so,Si) = SoSi[fi-2\{so, Sl) , 

(B.19) 

[L“](to,ti) = ti[xi-i](to,ti) • 


Proof. The integrated Legendre polynomials Li{si) for t > 2 vanish at si = 0 and si = 1 (see (2.17)). 
Hence, they must take the form 

Li{si) = Sl(l - Sl)V’i- 2 (si) , 

for some polynomial 'fi -2 G After homogenization one obtains the desired result. 


[L^] ( 50 , <si) — Li 


Si 


(so + Si)* — 


Si 


So + Si 7 So + Si 

SoSl / ,2 / f 


1 - 


Si 


fi-2 


Si 


So + Si 7 ^So + Si 

(so + Sl)^V'i_2(--)(so + Sl)*“^ = SoSl[V'j_2](so,Sl) 

Vsn + Si 7 


(so + Si)* 


where ['fi- 2 ] £ "P* ^(so,si) is the homogenization of V’i- 2 - The result involving L"(ti) follows using 
exactly the same reasoning and that it vanishes at ti = 0 (see (2.26)). □ 

"This result is not limited to Legendre and Jacobi polynomials and, as reflected in the proof, applies as well to the general 
polynomial families presented in Appendix A (see (A.2)). 
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Remark. Let Li and -L“ be the integrated Legendre and Jacobi polynomials of order i >2 and j > 1. Then 
there exists a homogeneous polynomial Xj-i] ^ si, S 2 ) such that 


(so, Sl,S2) = [Li]{so, Sl)[L"](so + Sl,S2) = SoSlS2[V'i-2, Xi-l](sO, Sl, ^ 2 ). 
Remark. Let Lj^ be the integrated Legendre polynomials of order k > 2. Then, 

[^fc] (m ([V’fc-2] ([V'fc-2] ( 1 , ^) G ^ 

where [^^^- 2 ] G si) is the homogeneous polynomial in Lemma 5. 


(B.20) 


(B.21) 


B.l Shape Functions 

For ease of referenee, note the deformed underlying spaee for m = 0 is 

VW-P = {</.€ QP/^P : V</.G X QP/-^’P-^ x 


B.1.1 Vertices 


The shape funetions for the vertiees are given by (9.11), and they are preeisely the 3D pyramid affine 
eoordinates. Take for example any quadrilateral vertex, say vi, and the top vertex ^ 5 , whieh merits its own 
attention. Their eorresponding shape funetions already satisfy the traee properties as diseussed when they 
were defined, so fhey lie in fhe lowesf order eompafibilify spaee Also, fhe verfex funelion for vi 

salisfies 

(F{x, y, z) = Ai(x, y, z) = ^ gM.i ^ 


, -(i-y) 

/ 1 + 2 

^4>'{x,y,z) = V\i{x,y,z) = [ 

^ ( 1+^)2 


■ S+’°' 


e 2 ; 


1,0.0 


(B.22) 




The same holds analogously for fhe ofher fhree quadrilaferal verfiees. Similarly, fhe lop verfex funelion 
salisfies 


<Fix, y, z) = X^ix, y,z) = , 

V<j/{x, y, z) = VA 5 (x, y, z) = 


0 

0 

1 

(i+Jp 


0 

0 

r)0,0,0 

^2 


^ 2 ( 


gO.l.O. 

1 , 0,0 


(B.23) 


Q+'“. 

;(0),1 


Therefore, all deformed verfex funetions lie in fhe lowesf order underlying spaee Voo'^’ , so lhal fhe original 
shape funetions are in fhe lowesf order spaee 


B.1.2 Edges 

Mixed Edges. The shape funelions for mixed edges are presented in (9.12). For every mixed edge 
fhey are labeled as fl, where i = 2,... ,p. As diseussed when fhey were defined, fhey salisfy fhe desired 
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trace properties, so they lie in the compatibility space To see they also lie in the underlying space take 

as an example mixed edge 12, and note that 




2,2,2 

2 


A{x,y)£Q'^''^{x,y) 


V(j)l{x,y,z) = 


(1 + 2 )' 


VA{x,y) +A{x,y) V 


(B.24) 


£Qi-l,i,OxQi,i-l,Ox{ 0 } 


(1 + 2 )' 
e{o}x{o}xQ°;°i° 


G s: 




S' 


i + l 


Here it was used that [Li] is a homogenized polynomial, so in particular (2.32) holds. Therefore, the de¬ 
formed shape functions are in so that the shape functions lie in Naturally, the same calculations 

hold for all other mixed edges. 


Triangle Edges. The shape functions for triangle edges are defined in (9.13), and are labeled as 0?, 
where i = 2,... ,p. They satisfy the desired trace properties, so they lie in the compatibility space T^^^’*. 
As an example take triangle edge 15. Then it follows 


(t>t{x,y,z) = (l)f{Xi5{x,y,z)) = [L,i\ A i+l T^) = -(T^ [Li]{{l-x){l-y),z) eQ:. 


i 2 , 2,2 

i 


A{x,y,z)GQ 




q: 


•i——1 . 


V(j)l{x,y,z) = 


( 1 + 2 )' 


VA{x,y,z) 


+A{x,y,z) V 


(B.25) 


ggi-2,i-l,i-lxQi-l.i-2,i-lxQi-l,i-l.i-2 


( 1 + 2 )' 
e{ 0 }x{ 0 }xQ°;“i° 


G Sr 


i + l 


where [Tj]((l — x)(l — y),z) = A{x,y,z) G due to Lemma 5, and where it is used that 

C Hence, the shape functions are in the correct underlying space and they 

belong to Analogous calculations hold for the other triangle edges. 


B.1.3 Faces 


Quadrilateral Face. The quadrilateral face functions are defined in (9.14) and identified as (plj, for 
i = 2,... ,p and j = 2,... ,p. Again, fhe funcfions are known fo safisfy fhe desired frace properfies, so fhey 
lie in fhe compafibilify space T^^)’®®, where n = max{i, j}. To see fhey are in fhe underlying space, nofe 

(l)lj{x,y,z) = po{j^)4>?j{m{x),yoiiy)) = j^(j)fj{m{x),fioi{y))eQ\^’° c Q", 

'-V-' 

A(a;,j/)eS'’l(a;,j/) 


,n—l.n.n —1 ■ 


(B.26) 


,n,n —l,n —1 


q: 


n,n,n — l 
n + i 


'^(l>ijix,y,z) = VA{x,y) +A{x,y) V jA ^ 

eS'-i.2.0xQi,2-i.Ox{o} e{o}x{o}xQ°-“’° 

where fhe inclusions of fhe fype ^ C are used repeafedly. Nofe fhaf fhe facfor 

PoiiA) required in order for fhe funcfion fo be in fhe correcf space (see fhe firsf fwo componenfs of 
fhe gradienf). If follows fhe shape funcfions belong fo where n = max{z,y}. 
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Triangle Faces. The triangle faee funetions are defined in (9.15) and are labeled as cjilj, for i > 2, 
j > 1, and n = i + j = 3,... ,p. They satisfy the traee properties, so they lie in Take for instance 

face 125, and observe the functions also lie in the corresponding underlying space, because 


V4j(x,y,2:) 


n — 1, n — 1, n — 2 


yo{y)(t)ij{i'oi 2 {x, z)) = yo{y)[Li, Lj]{l - x,x,z)eQn’ 

VA{x,y,z) 


A{x,y,z)&Q' 


(l+z)" 

-- 

^Qn —2,n — 1,71 —2y^Q'n—l ,n —2,n — 2y^Qn — l ,71 —1 ,rt —3 


+A{x,y,z) V 


(1+^) 
e{0}x{0}xQ°'“f 


G S 


.n—l,n,n —1 ■ 
■n 

in,n —l,n —1 


(B.27) 


s: 


n + l 


Here, [Lj,L"](l — x,x,z) = A{x,y,z) G Q" ^ by (B.20). Hence, (plj is in the correct space 

^(0)'"^, where n = i + j. The same follows for the other triangle faces. 


B.1,4 Interior Bubbles 


The interior bubbles are defined in (9.16), and identified as 0^^^, where i = 2,... ,p, j = 2,... ,p and 
k = 2,... ,p. They satisfy the vanishing properties along the whole boundary so they are trivially in 
for n = max{i, j, k}. They also belong to the underlying spaces since 


4>\j{x,y,z) 

'^4>\j{x,y,z) 


(l>k (mCily)) (i)ij{m{x),poi{y)) G 








A{x,y)&Q*Z[x,y) 


Biz) 


VA{x,y) +A{x,y) VBiz) G 
e{o}x{o}xQ°-°f-i 


n —l.n.n—1 ■ 


n.n—l,n—1 


eQi-l,i,OxQij-l,Ox{0} 


Q: 


n,n,n —1 
n+l 


(B.28) 


Here, (B.21) was explicitly used. It follows all the interior bubbles lie in where n = max{i, j, k}. 


B.2 iT(curl) Shape Functions 


Recall the deformed underlying space for m = 1 is 



{EeQ 


p—l,p,p 

P+1 


X Q 


p,p—l^p 

p+1 


X Q 


p,p,p—l 

p+1 


VxEgQ 


p,p—l,p—l 

P+2 


X Q 


p—l,p,p—l 

p+2 


X Q 


p—l,p—l,p—l 

p+1 


}■ 


B.2.1 H (curl) Edges 

/^(curl) Mixed Edges. The mixed edge shape functions are given in (9.17), and are labeled as Ef for 
i = 0,... ,p — 1. As discussed before, the trace properties are satisfied, meaning that they belong to 
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To see they are in the underlying space, consider edge 12 and note that 


E!{x,y,z) = no{y)E,^ 


E/ 1—x 


(1+^) 


soW 5rj(a;) 


i+2 I 


gQ 0 , 0 , 0 xg 0 , 0 , 0 x{ 0 } 


= (i+l )»+2 Mo(j/)[-Pi]('So(a;),si(a;))(go(x)Vsi(a;)-gi(x)Vso(x)) £| 


G{0}x{0}xQ^’^’0 




/ 

^i + l,i,i 
' ^i+Z 

xA{x,y)e[ 

^i,i+l,i 

^i+Z 

- \ 



i,i+l,i+l , 
i+2 

i + l,i,i+l 
i+2 

i + l,i + l,i 


' 2i+2 


(B.29) 


Hence, the pullback of the shape functions lies in so the original shape functions lie in as 

desired. The same calculations hold for the other mixed edges. 


iT(curl) Triangle Edges. The triangle edge functions are defined in (9.18), and are identified as E|, 
where f = 0,... ,p — 1. Again, fhe frace properfies are satisfied, so fhaf fhey are inside Take for 

insfance edge 15. Making use of (4.9), (B.22) and (B.23), notice fhaf 


Sl(^) 

Et{x,y,z) = Ef{\i{x,y,z),\;^{z)) = j^^^^Ef{{l - x){l - y),'^) 

ggO.i.ixQi.o.ixQi.i.o 


■a 


i^i+l,i+\ . 
i+2 


= (i+lf+2 [-Pi](go(a;, y), si( 5 ;))(so(a:, 2/)Vsi(2:)-si(2:)Vso(a:, y)) G S^a 


i + l,i,i+\ 


■Q 


i + l,i + l,i 
i+2 


^q0,i,0xq1,0,0xq1,i,0 e{o}x{0}xQO’0’'' 


VxE- (x, y, z) = [Pi] (so(x, y), si(z)) ^ VXi{x,y,z) x VA5(z) G | 

eQ3’°’°xQ“’^’°x{0} 


i + 3 


i+2 


(B.30) 


Therefore, fhe shape functions belong fo and as a resulf lie in The same reasoning is 

attached fo fhe ofher friangle edges. 


B.2.2 H (curl) Faces 

H (curl) Quadrilateral Face. The fwo closely related families of quadrilaferal face funclions are pre¬ 
sented in (9.19) and (9.20). They are labeled as Efj, for i = 0,... ,p — 1 and j = 2,... ,p. As usual, 
fhe funclions are shown fo salisfy fhe desired frace properfies, so fhey lie in fhe compafibilify space 
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where n = max{i + 1, j}. They lie in the underlying space because 


E: 





n+1 


A(x,j/)eS®'J’°xQt.i,Ox{o} 

e{0}x{0}xQ“’°'° 

VxEl^{x,y,z) = j^^ y X A{x,y)+ xA{x,y)e[ Q 

' ' ' ^ V S: 


-«n,n —l.n—1 


(B.31) 


e{ 0 }x{ 0 }xQ"-i."-i.o 


eQj‘'°xQ‘'^''“x{0} 



where the inclusions of the type ^ Qn+i~^ — 2n+i repeatedly. Notice that the factor 

/io(x^)^ was required in order for the function to be in the correct space (see the last component of the 
curl). It follows the shape functions belong to where n = max{i + 1, j}. The calculations are 

invariant to permutations, so the result holds for both families. 


H (curl) Triangle Faces. The two closely related families of triangle face functions are defined in (9.21) 
and (9.22). They are identified as Ef^, for ^ ^ 0, j > 1, and n = i-\- j = 1,... ,p — 1. They satisfy the trace 
properties, so they lie in Take for instance face 125. Using Lemma 5 it follows that 


Eij{x,y,z) = iJ,o{y)Efj{doi 2 {x,z)) = po{y)[Pi, L‘j]{doi 2 {x, z))EQ{iJoi{x, z)) 


_ ^n — l,n — 1 ,n — 1 
^^n-1 




eS: 


n,n,n^Qn,n,n^^n+l,n+l,n—1 


n+2 


'n+2 


xQ: 


VxElj{x,y,z) G Vx 


, q: 


Q n,n,n 
n+2 

Q n,n,n 
n+2 

n+l,n+l,n —1 
n+2 


Q s: 


n+2 
n + l,n,n , 

n+3 

n,n+l,n 

n+3 

Q n,n,n 
n + 2 


■q: 


= [Pi,Xj-i]ii^0i2ix, z)) fio{y)i^2{x, z)EQ{i 7 oi{x, z)) £ s: 


S, 


n,n+1,n+1 , 
n + 2 

n + 1,n,n+1 
n + 2 

n + 1,n + 1,n 
n + 2 


(B.32) 


Note the calculations hold regardless of permutations in i'qi2{x, z), meaning that both families lie in the 
underlying space Hence, the shape functions Elj are in the correct space where n = i+j. 

The same follows for the other triangle faces. 


B.2.3 H (curl) Interior Bubbles 

Family I: The first family is defined by (9.23) and it is the gradient of the bubbles. They are 
labeled as where i = 2,... ,p, j = 2,... ,p, and k = 2, ... ,p. They satisfy the trace properties and 
lie in for n = max{f, j, k}. As gradients of they lie in the underlying space. Indeed, from 
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(B.28) it follows 


Efjkix,y,z) = v(^(p^{lloi{j^))4>?j{mix),i2oi{y))"j el a 


-»n —l,n,n —1 
iln 

-jn,n —l,n —1 
/-«n,n,n —1 


VxEfji^{x,y,z) = 0e QZ+l’"-’" ^ 


-«n —l.n —l.n —1 


Hence, the interior bubbles lie in for n = max{z, j, k}. 


(— I — 

^ 1 ^n+l 

V 


(B.33) 


Families II and III: These families are presented in (9.24) and (9.25), and are identified as with 
i = 0,... ,p — 1, j = 2,... ,p, and k = 2,... ,p. They satisfy the vanishing trace properties and belong to 
r(i)>"^, for n = max{i + l,j, k}. Using (B.31) and (B.21) it follows 


-»n —l,n,n —1 


--V- ^ -V -^ V 0 / V / 

^ -^n+i 

tz.f nls//^ .34^ 

G|0}X{0}XC2^ 2 , ^n,n-l,n-l ^ 

/ ^n + 2 \ 

'^xE'°jk{x,y,z) = B{z) V X A{x,y) + VB{z) x A{x,y) e [ • 

^ ^ \ ^n-l,n-l,ri-l / 

6 {o}x{o}xQ"-i>»-i>o eQ" 4 ’ 2 *=-Us^’:^’*-U{o} ^"+i 

The analysis is independent of permutations of the entries, so it holds for both families. Therefore, the 
functions of both families are elements for n = max{z + l,j, k}. 


^n + 2 

-«n —l.n.n —1 


-»n — 1 , n — 1 , n — 1 


Family IV: The family is shown in (9.26), with the functions being identified as E^j, for i = 2, 
and j = 2,... ,p. They have vanishing trace, so they lie in for n = max{i, j}. Moreover, 


E'^jix^y^z) = (/>°(^oi(a;),Moi(y)) ^(iq^ 


/ 0 \ (QT+A^ 

G 0 J C 


A(x,y)(iQ-’i{x,y) B(^)e{0}x{0}xQ°’°'“ 


VxEij{x,y,z) = VA{x,y) xH(2:)g(c 

eQi-i,i.OxQi.i-i.Ox{o} ' 

Hence, the functions belong to for n = max{i, j}. 


/ \ 

, / 

Qir+f’° 


' “ \ 

Vo/ 



Q n,n —l,n —1 
n +2 

S n—l,n,n —1 
n +2 

.-«n — l.n—l.n —1 


(B.35) 


B.3 i?(div) Shape Functions 

Recall the deformed underlying space for m = 2 is 
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B.3.1 iJ(div) Faces 


Ff(div) Quadrilateral Face. The quadrilateral faee funetions are defined in (9.27), and they are labeled 
as for z = 0,..., p — 1 and j = 0,...,p—1. The functions satisfy the desired trace properties, so they 
lie in the compatibility space where n = max{z +1, j +1}. They lie in the underlying space because 


ylj {x, y, z) = (^01 (x) , m (y ))' 


A(a;,j;)g{0}x{0}xQ®>J>0 


Q n—l,n —1,0 
3 


c s: 


-»n,n—l,n—1 , 
^n+2 

-»n —l,n,n—1 
^n + 2 

-.n —l,n —l,n 


V-V^{x,y,z)= • Aix,y) G C 


n—1,12—1,71—1 


(B.36) 


E{0}x{0}xQ°’°'“ 


Notice that the factor was required in order for the function to be in the correct space (see the 

divergence). It follows the shape functions belong to where n = max{z + 1, j + 1}. 


H (div) Triangle Faces. This is by far the most difficult construction (and proof). Since only a sketch 
was shown before, more details will be given here in relation to those functions. As an example consider 
face 125. 

The typical approach is to use the ancillary functions with appropriately chosen affine coordinates as 
their input, and possibly some blending factor to ensure the trace properties are satisfied. Indeed, this 
directly yields all the shape functions up to this point. However, the two most sensible options for the shape 
functions fail to be even in the lowest order space. 


MvWooi^oMx, z)) 

'y ■ [yo{y)yooim2ix, z))'^ 

jI^yooi^i25{x,y,z)) 



(B.37) 


One option is to depart from the spaces defined by Nigam and Phillips (2012), by adding elements to 
^(2),p However, this has the negative effect that the discrete space, also has to be modified to 

satisfy the exact sequence property. Thus, one could have scenarios where the lowest order discrete space 
for is the constants plus other functions, as opposed to simply the constants. Although valid, this deviates 
from the other elements, where the lowest order discrete space is always the constants. Moreover, this 
option adds degrees of freedom to the construction. Hence, modifying the spaces and is a last 

resort, and should be avoided when possible. 
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A second option, is not to use the ancillary functions. However, using the ancillary functions directly is 
not a capricious decision. The main reason is that they automatically guarantee that no “illegal” derivatives 
are present in the divergence of the shape functions. Indeed, the traces should be of the form of 2D 
functions for a given face, so intuitively, elements of should be involved in the expressions for the shape 
functions. However, the divergence of the shape functions cannot involve derivatives of those traces, since 
they are elements of L^, which in general do not have derivatives. Indeed, the traces for the triangle faces 
in our shape functions involve combinations of Legendre and Jacobi polynomials. Pi and which are the 
representatives of L?'. As expected, up to now the divergence of i7(div) face functions has not involved 
any derivatives of P* and P". This is not a triviality, since in general one would expect the derivatives to 
be there, but due to Lemma 4, the derivatives disappear if the ancillary function 14^ is utilized. Therefore, 
amongst other reasons, it is highly desirable to use the ancillary functions. 

With these facts in mind, the message is then to at least try to persist in the use of the ancillary functions 
without modifying the space. The key to making this possible relies in observing closely the lowest order 
space. Here, the first order shape functions are known explicitly for each face and have been deduced and 
presented by Gradinaru and Hiptmair (1999) and Nigam and Phillips (2012) amongst others. For face 125, 
the pullback of that shape function is 




1 

(T+iF 


0 

-(i-y) 

I- 


■ ■ 
qO.i.o 


V-VoQ{x,y,z) = G Q 4 


0 , 0,0 


(B.38) 


As a result, any proposed shape functions should match this expression in the lowest order case. Even 
though this was not the case for the previous attempts shown in (B.37), a close observation reveals that a 
linear combination of those attempts does lead to the desired lowest order shape function. That is,'^ 


1 

2 



(Foi2(x,z)) 


^^00 (^125(2:, y,z))) 


2(1+F" ( ^ Q + 2(TTiF ( ) 


(B.39) 


Immediately, this suggests the higher order expression for the shape functions. 


Vij{x,y, z) = l(^Ho{y)V^{doi 2 {x, z)) + -^V^{Xi25ix,y, z))^ , 

V-V^j{x,y,z) = l(yfio{y) ■ V^{i7ou{x, z)) + -^V ■V^{Xi 25 {x,y, z)) (B.40) 

- J^'^yoiy) ■V^{Xi 25 {x,y,z))'^ , 


where f > 0, j > 0 and n = i + j = 0,... ,p — 1. Here, it was used that V • Vj^{i'oi 2 ix, z)) = 0 by (7.14). 
Clearly, no derivatives of the representatives are present in the expression for the divergence, since they 
are not present neither in the terms involving directly, nor in the terms involving V • in view of 

'^In fact, the case z), X 2 {x, y, z), As( 2 )) + Vqo (Ai(a:, y, z), ui{x, z), As( 2 ))) also gives the desired lowest order 

shape function. However, the higher order version of this expression unfortunately is not in the space . 
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f (2) n+1 

Lemma 4. It remains to show the traee properties, and that the are in the underlying spaee V(» , 

where n = i + j. 

In faet, eaeh of the two eomponents of the shape funetion satisfies the traee properties. To see this, note 
that Vj^{i 7 oi 2 {x, z)) has entries independent of y, meaning that only the seeond eomponent is nonzero. This 
eomponent is tangent to faees 235, 145 and 1234, so that the normal eomponent vanishes on those three 
faees. Meanwhile in the opposite faee 435, yo{y) vanishes, so the shape funetion also vanishes, while at 
the faee itself yoiv) = and the shape funetion (in the pyramid {C,y,C) eoordinates) takes the form of a 
2D triangle faee funetion, as desired. Therefore, the eomponent yo{y)Vj^{ uqi2{x, z)) satisfies the traee 
properties. For the other eomponent, notiee that 

y, z)) = [Pi, P^]{Xi 25 (,x, y, z)) (^^o(f/)^oo (^ 012 ( 2 :, z)) + X5{z)Vyo{y) x (Foi(x, z))^ . 


Henee, the seeond eomponent is further deeoupled into two terms, with the first one already satisfying the 
desired properties by the previous analysis. The seeond term vanishes at all faees as required, beeause 
V^o(y) is normal to faees 125 and 435, while z)) is normal to the faees 235 and 145, meaning 

that Vy,o{y) x Eq{Pqi{x, z)) is tangent to all the triangle faees, so its normal traee vanishes. Finally at the 
quadrilateral faee it vanishes beeause A 5 (z) = 0 there. It follows yo{y)~^V^^{Xi 2 ^{x,y,z)) satisfies the 
traee properties as well, and so G p( 2 ),n+i ^ ^ 

To prove the shape funetions lie in the underlying spaee, simply reeall that [Pi, P“] G si, § 2 ) 

is a homogeneous polynomial, so it suffiees to assume it is a monomial [Pi, P^]{so, si, S 2 ) = of 

order n = i + j = a + b + c. Caleulating explieitly for faee 125 gives. 


yij{x,y,z) 

'y-Vij{x,y,z) 


2 ( 1 + 5;)"+3 




\q 


n+l,n,n ^ 

n+3 

n,n+l,n 

n+3 

n,n,n + l 

n+3 


2:^(1 - x)“x^((n + 2 - 2:)(1 - 2/)“+^ + (1 + z)) 

2(1 + z)^+* 


(B.41) 


In the expression for the divergenee, the power of 2; is either c or c + 1. If a + 6 = 0, then c = n and 
V-V^j{x,y,z) = G Qn+f - ^ a + b > 1, then c + 1 < n and again V-V^j{x,y,z) G 

Therefore in all eases, V^j{x, y, z) G where n = i + j, and as a result G Naturally, 

an analogous result holds for all other triangle faees. 


This part eoneludes with the observation that the expressions in (B.40) have faetors 


m(j/) 


and 




whieh appear to be singularities (on faee 435). However, they are not real singularities, sinee 


hV/^iyso, fisi, S 2 ) = [Pi, P^]{yso, fisi, S 2 ){yVQ^{so, si, S 2 ) + S 2 V yx Eq{so, si)) , 

/ \ / A N (B.42) 

V- [j^Vifiyso, fisi, s2)j=[Pi, P^Kyso, fisi, S2)Vy(ii+j+3)E^{so, si)xVs2-'Foo (so, si, S2)j, 

for f > 0, j > 0 and n = i + j = 0, ...,p — 1. Here, in the expression for the divergenee the term 
//(Vso X Vsi) • Vs 2 is not shown, beeause it is assumed that so + si + •52 = 1, sinee this is the property F 012 
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satisfies. This expression can be very useful from a computational standpoint, because the benign singularity 
can be a problem when computing the shape function at interior points very close to face 125. Indeed, if 
orientations are not being taken into account, the expression above is sufficient to avoid those problematic 
terms. Nevertheless if orientations are considered, the solution is more technical. We present an approach 
that is convenient from a computational point of view. First define. 


k{o) 


0 ifo = 0,l,2, 
1 if 0 = 3,4, 5. 


(B.43) 


Then define the following funetion, 

'^ij (^0? ^1?'^2;^0? ] (^o (^oq ’ ^ 2))) (*^oo (^0’*^1 ’ ^2))) 

V-Vif{so, si, S2, F, 00 , o) = [Pi, P“](o-^(cj^(/rso, S 2 )))V/r (B.44) 

• ((f+j+3).E^(cj^(„)(cj^(„^)(so,si)))xVs2 

for i > 0, j > 0 and n = i + j = 0,... ,p — 1. Then, the orientation embedded shape functions for each 
face (see (9.28)) are 

0 = 2 0))) + ^^^(^^ 012 ( 0 ,0, Fc(i^), oq, o)) , 

; .. X (B.45) 

V-F4(C,0 = ^VMOi^)-^.fOo^«fe(ea,0)))+V-y,f(Foi2(0,0,/Uc(i^),oo,o)), 

where i > 0, j > 0, n = i + j = 0,... ,p — 1, {a, b) = (1, 2), (2,1), c = 0,1, and where oq is chosen 
such that o'^(Foi 2(0) 0) is the locally oriented triplet representing that face. This is necessary because S 2 
is treated differently than (sq, si) in the definition of 14^- 


B.3.2 H (div) Interior Bubbles 


Families I and II: These families are presented in (9.29) and (9.30) and are the curl of H (curl) interior 
bubbles. They are identified as 14^^, with i = 0,...,p—l,j = 2,...,p, and k = 2,... ,p. They satisfy the 
vanishing trace properties and belong to for n = max{i + l,j, k}. Using (B.34) it follows 


y^jk{x,y,z) = Vx(^j^(j)f{iloi{j^))E^j{fIoi{x),fIoiiy))^ 


^ ^n,n — l,n — l^n — l,n,n — l^n — l,n — l,n — l 
^^n+2 ^^n+2 ^^n + 1 


V-V^j{x,y,z) 


0 £ Q 


n+3 


G 


Q n,n—l,n—1 ^ 
n + 2 

Q n — l,n,n—l 
n + 2 

Q n — l,n — l,n 
n + 2 


(B.46) 


The calculations are invariant to permutations of the entries, so it holds for both families. Therefore, the 
functions of both families are elements where n = max{f + l,j, k}. 
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Family III: This family is defined in (9.31) and it is the curl of H (curl) interior bubbles. The functions 
are labeled as V^, for i = 2,... ,p, and j = 2,... ,p. They have vanishing trace, so they lie in for 

n = max{i, j}. Moreover, using (B.35) gives 


y^jix,y,z) = VX fIoi{y))V 

Hence, the functions are in for n = max{t, j}. 


n+l 
1 ^- 1,.?',0 
■n + l 


n,n —l,n —1 , 



(B.47) 


Family IV: This family is shown in (9.32) and the functions are identified as 17^^, for i = 0,..., p — 1, 
j = 0, ... ,p — 1, and k = 2, ... ,p. They satisfy the vanishing trace properties and lie in the compatibility 
space where n = max{t + 1, j' + 1, k}. Using (B.36) and (B.21) it follows 


yijki^,y,z) 

V-Ui5fc(x,|/,z) 


= (Moi(if^)) Ui°(/roi(x), /2oi(y)) e 


Q. 


0 ' 
0 

n —l,n —l,n—1 
n+2 . 


c 


A(x,U6{0}x{0}xQ*.J.o 

i,j,k—l f- ^n—l,n—l,n—1 


= vm ■ Aix,y)€Q]:ircQ:-;, 

e{0}x{0}xQ°^’''-i 


Q 

Q 

Q 



(B.48) 


Therefore, the functions belong to for n = max{i + l,j + 1, k}. 


Family V: This family is defined in (9.34), with the functions being labeled as V^, for i = 2,... ,p, 
and j = 2,... ,p. They vanish at the pyramid boundary and belong to the compatibility space where 

n = max{i, j}. Using the expression in (9.33), it follows 


■q: 


V!j{x,y,z) = {^)^ 1 VA{flQ^{x),floi{y),po{i^)) G s: 


^l + Z' 


r< 

' ij 


n,n—l,n—1 , 
n + 2 

n —l.n.n—1 


l+z^ 


V-V^^ix,y,z)= V(^) 


n-l . A(.r. ^ n*-14-T^-2 


1 + 2^ 


. q: 


^n + 2 

n —l,n —l,n 


n + 2 


(B.49) 


^ l,n—l,n—1 

- ^n+3 


e{0}x{0}xQ°’°'"-2 


As a result, the functions are in for n = max{z, j}. 


Families VI and VII: These families are presented in (9.36) and (9.37), with the functions being 
labeled as and respectively, where i = 2,... ,p, and j = 2,... ,p. They satisfy the trace properties 
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and therefore are in the spaee and respeetively. Using the expression in (9.35), it follows 


Vi'ix, y, z) = {T^y ^ U.-(/Ioi(x), /ri(y), ^o(t4t)) G 


a+2 


ggOAi-i A(x,y)GQ^°’°x{0}xQi-^’°’° 


1+z' 


^i+2 

^i+2 

i+2 


V-Vy{x,y,z)= V(^) 


2—1 


A{x,y) G C 2,^3 


2 — 1 , 2 — 1 , 2—1 


e { 0 } x { 0 } xQ “’°''-2 


(B.50) 


Symmetrie arguments apply to the seventh family Vj’{x, y, z), meaning that the funetions are in and 
^(2)J respeetively. 


B.4 Shape Functions 

Reeall the deformed underlying spaee for m = 3 is 

v(3),p _ qP-1,P-1,P-1 

>^00 — ^ p +3 


B.4.1 Interior 

The interior shape funetions are defined in (9.38) and identified as forz = 0,...,p — l,j = 0,...,p—1 
and k = 0,... ,p — 1. They trivially satisfy the eompatibility properties, sinee there are none to satisfy, 
meaning they lie in for n = max{f + l,j + l,k + 1}. In this ease the faetor making the expression 

eoordinate free is {Viyi{x, z)xVi/i{y, z)) -V It follows, 

i’'ijk{x,y,z) = [Pi\{m{x))[Pj]{poi{y))[Pk]{m{T^)) g 21+4 ^ • (^-51) 

'-V-' 

Henee, the interior funetions lie in for n = max{f + l,j-|-l,/c + l}. 
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C Integration 

C.l Coordinate Changes 



Physical Element 



Figure C.l; Transformation from the master element to the physical element. 

At its core, the finite element method advocates carrying out integration over a master element domain 
instead of the original physical element. It makes the method very feasible from a computational standpoint. 
This involves a change of variables r : n —)• from the master element domain Q to the physical domain 
O, which is assumed to be known. This is illustrated in Figure C.l. Note the change of variables r is in 
general a nonlinear mapping. 

Indeed, consider a “physical” integrand / which is a function of variables in the different energy spaces 
and their differential form. These variables are in the physical system of coordinates. For instance, take 
E^, and to represent variables in Ff(curl), Ff(div) and respectively. Their corresponding 
differentials are x and However, it is their pullbacks to the master element domain, 

denoted with the subscript H, which are known, since the shape functions are defined in the master element 
domain.Making use of the appropriate pullback mapping for each of the variables as written in (B.8), this 
yields 


= jj ° (yM °EE^oT,iV^x E^) o r, V~ o r, (V^ • V~ ) o r, V’o o ciet( J^)dO 

detljr)^^^' d^tiX)V’tt)det(Jr)dn , 


(C.l) 


where is the Jacobian matrix of the transformation r : fl —)• fl. 

'^Unless the affine coordinates and their gradient are written in the physical system of coordinates, in which case one can simply 
substitute them in the expressions for the shape functions. This is due to the coordinate free nature of the shape functions. 

'"'in 2D, Vf 2 X is in L^, so the correct expression in the last line would be ) Vn x En instead of the 3D expression 

det(j ) -Jt Vn X En- In ID, En and Vn do not even exist, so they would be ignored throughout. 
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Now, the integration is at least in a well known master element domain. However, this is still a “difficult” 
domain over which to integrate (with the exception of the ID segment, the 2D quadrilateral and the 3D 
hexahedron). Hence, it is desirable to make one further change of coordinates to a “nicer” integration 
domain. This is denoted by rg : Q —)• D, where Q = (0,1) in ID, Q = (0,1)^ in 2D and Q = (0,1)^ 
in 3D. The transformations for each element are nicely depicted in Figure C.2. Some readers may have a 
strong preference for the integration domains Q = (—1,1) in ID, Q = (—1,1)^ in 2D and Q = (—1,1)^ 
in 3D. If that is the case, then simply make the substitutions x = y = and 2 = in those 
expressions shown in Figure C.2. 


^ = x 

v = y 


- 

Quadrilateral 


2: = ^ 

y = v 


y 


--^ X 

Integration Square 


5 = 2(1 - y) 

y = y 




Triangle 
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The original integral in (C.l) finally beeomes 


f{<Pii°TQ,{JT'^'^n4>n)oTQ,{J^'^En)oTQ,{-g^^JrVnxEn)oTQ, 

■^Q ^ (C.2) 

° ° (ditpo'/’tt) o rQ)det(J,Q) det(J,)dQ , 

where Jtq is the Jaeobian matrix of the transformation tq : Q ^ 

C.2 Fast integration 

To aetually ealculate the integral, the typical approach is to use Gaussian quadrature. In Q (or Q) the 
quadrature points and weights are well known and taken from the literature, and this is part of the reason 
why integration over the physical domain was reduced to integration over Q in (C.2). 

However, as the number of spatial dimensions N increases from ID to 3D, the cost grows quickly with 
p. Indeed, to construct a typical finite element stiffness matrix, integrals usually reduce to the form 



where / = 0{p^) and J = 1,..., 0{p^). The cost to integrate each term is 0{p^) as well, because 

there are p + 1 quadrature points in each spatial dimension. Hence, with a straightforward implementation, 
the cost to integrate all terms is 0{p^^), so it is 0{p^) in ID, 0{p^) in 2D and 0{p^) in 3D. This constitutes 
a problem for > 2 and high p, so it is highly desirable to improve the integration cost. 

Fortunately, the integration cost can be reduced if there exists a decoupling of either uj or vj in x, y 
and z (the variables after transforming to Q, not necessarily the variables of the physical or master element 
domain). Assume the decoupling is in vj, where it takes the form vj{x,y) = nf {x,y)v^ (y) in 2D and 
vj{x, y, z) = {x, y, z)v^-^ (y, 2 ;)n|^ {z) in 3D, with jx,jy, jz = 1, • • •, 0{p). Then, by reorganizing the 

operations and storing some coefficients, the cost is reduced to 0(p^) in 2D, and to 0{p'^) in 3D. Some 
of the details are in Demkowicz et al. (2007). With the shape functions presented in this text, regardless 
of the element shape and the the associated topological entity, such a decoupling is to be expected, so this 
acceleration to 0{p^^~^^) is possible. This technique based on a tensor product decoupling is typically 
cdlXeAfast quadrature. 

Naturally, there are other fast integration techniques different from the fast quadrature described above 
which might also be applicable, but further research is required. 
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D Verification 


One of the most important tests is to numerieally eonfirm the polynomial approximability properties of the 
spaees spanned by the shape funetions. Coupled with the exaet sequenee property of the diserete spaees, 
this ensures all well known interpolation inequalities. More speeifieally, for an affinely transformed master 
element mesh, let be the span of the basis funetions of order p (being eomposed pieeewise by shape 
funetions), and similarly with Q^, and for the spaees H{div) and respeetively. Then, 

one has to eheek that VP C WP, (pP-^f C Qp, (pP-^)^ C VP and pP-^ C YP. 

To do this, first eonsider an arbitrary rt in a given energy spaee U approximated by a diserete spaee Uh- 
Clearly, 

u G Uh dist{u,Uh) = min \\u — Uh\\u = d . (D-1) 

Uh&Uh 

Henee, given u G U the task is to eompute dist(tt, Uh) = ||rt — Uh\\u^ with Uh G Uh being the element 
where the minimum is attained. Fortunately Uh is eomputed from a variational problem equivalent to the 
projeetion (distanee) problem. It is, 

I Uh ^ 1 

\\u - UhWu = dist {u,Uh) < (D.2) 

yb{uh, Vh) = {uh, Vh)u = {u, Vh)u = ^u{vh) yvh G Uh- 

Naturally, the inner produet is different depending on the energy spaee U. They are. 


{4>i, 4>2)m 

{Ei,E2)h{ curl) 
{Vl,V2)H(d iv) 

(V'l,V’2)l,2 


[ (</>l(/>2 + V(/)1 • V02)dn , 

Jn 

[ {El ■E2 + {Vx El) ■ (V X E2))dn , 
Jn 

[ (Ci-V2 + (V-Ci)(V-V2))dF!, 

Jn 

[ 'ilJi'ip2dO,. 


(D.3) 


Then, the task is to determine whether eaeh element tt of a monomial basis for the polynomial spaees 
in question lies in Uh- This is aehieved by solving the variational problem in (D.2) and eheeking that the 
relative error range of maehine zero. Thus, for example to ensure that 

pp C WP one must eheek that all monomials of the form for i + j + k < p lie in WP, the span of the 

basis funetions. Similar proeedures hold for H (curl), H (div) and L^. These tests are ealled polynomial 
reproducibility tests. 

The polynomial reprodueibility tests are sueeessful when using the eode assoeiated with this work. The 
tests are done on a series of meshes, ineluding a four element hybrid mesh with one element of eaeh type, 
as depleted in Figure D.l. By doing this on a mesh, there is the additional value of implieitly verifying 
eompatibility of the shape funetions aeross the boundaries of the elements. Indeed, it should be eheeked that 
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Figure D.l; A hybrid mesh used to verify polynomial reproducibility. 

the polynomial reproducibililty tests pass under all possible orientations of each face and edge in the mesh, 
which is the case for our code. 

Another convenient test is to verify some aspects of the exact sequence property of the discrete spaces. 
For this, consider a fixed element and the discrete spaces W^, Qp, and conforming to H^, 

H{div) and respectively. The discrete spaces are precisely the span of the corresponding shape functions. 

Hence, for example consider an shape function cf) G W^. Then the idea is to confirm numerically 
fhaf V(/> G QP. This is done as described in fhe polynomial reproducibilify fesfs, where fhe compufed 
projecfion of V(/> fo Qp is E^, which is given as a linear combinafion of fhe Tf(curl) shape funcfions 
E G Q^. Therefore, one should obfain fhaf fhaf is in fhe range of machine zero. Moreover, 

^ l|V0||H(curl) ^ 

one can addifionally check fhaf fhe coefficienfs of fhe linear combinafion for E^ make sense. For insfance, 
a (j) £ is originally an inferior bubble, fhen V4>is also an H (curl) inferior bubble and as a resulf is 
in fhe span of fhe H (curl) inferior shape funcfions (meaning fhe coefficienfs associafed fo H (curl) edge and 
face shape funcfions are zero). Similarly, ii(j) £ is an face shape function, fhen Vi?i) is in fhe span of 
fhe fhe Ff(curl) face funcfions associafed fo fhe same face and fhe iF(curl) inferior bubbles. Nafurally fhis 
applies fo ofher fopological enfifies and fo fhe differenf energy spaces. 

These verifications are successful when using fhe code fhaf supplemenfs fhis fexf. 
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E Tables 


E.l Polynomials 

Polynomials 

Legendre 

Shifted and scaled Legendre polynomials for t > 0 and x G [0, f]: 

Poix;t) = 1 
Pi{x; t) = 2x — t 

iPi{x; t) = (2i — l)(2x — f)-Pj_i(x; t) — {i — l)f^Pj_ 2 (x; t) for i > 2 

Integrated Legendre polynomials: 

2(2i — l)Lj(x; t) = Pi{x; t) — t^Pi- 2 {x; t) for i > 2 
Scaling differential of integrated Legendre polynomials: 

PLi+i{x]t) = Ri{x]t) = -l{Pi{x;t) + tPi-i{x;t)) for r > 1 
Homogenized polynomials: 

[Pi]{so,si) = Pi{si;so + Si) for r > 0 
[Ri]{so,si) = Ri{si-,so +si) for i > 1 
[Li](so,si) = Li(si; So + si) for i > 2 
V[Li](so,si) = [Pi_i](so,si)Vsi + [i?j_i](so,si)V(so + si) for i > 2 
where so and si are functions of some spatial variable in so Vsq, Vsi G for N = 1,2, 3. 

continued on next page 
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continued from previous page 

Polynomials 

Jacobi 

Shifted and scaled Jacobi polynomials for a > — 1, f > 0 and x G [0, f]: 

Po{x]t) = 1 

Pi{x; t) = 2x — t + ax 

ajP^{x] t) = bj{cj{2x — t) + a^t)P^_i{x] t) — djt^P^_ 2 {x] t) for j >2 
where 

aj = 2j'(j + a) (2j + a - 2) 

bj = 2j + a — 1 

^3 ~ + a — 2) 

dj = 2(j + a - 1) (j - 1) {2j + a) 

Integrated Jacobi polynomials: 

L"(x; t) = X 

L“(x; t) = ajP^{x, t) + bjtP^_i{x; t) — Cjt^P^_ 2 {x; t) for j >2 
where 

a . - _ i±^ _ 

“j (2j+a-l)(2j+a) 

=(2i+a-2)(2i+a) 

_ j—1 

Cj - (2j+a-2){2j+a-l) 

Scaling differential of integrated Jacobi polynomials: 

§^Lj+i{x]t) = Rf{x;t) = + fP“_i(x;f)) for j > 1 

Homogenized polynomials: 

[P"](fo, ti) = P"(fi; fo + h) for j > 0 
[Rj]{to,ti) = Rj{ti-to + ti) for j > 0 
[L'^]{to,ti) = L'^{ti;to + ti) for j > 1 
V[L“](fo,fi) = [P“_i](fo,ti)Vfi + [R^_,]{to,ti)V{to + ti) for j > 1 
where to and ti are functions of some spatial variable in so Vfo, Vfi G for N 




E.2 Ancillary Operators 


Ancillary Operators 

In all cases, sq, si, S 2 , to functions of some spatial variable in ] 



Operator 

Indiees 

Edge 

= [L*](so,si) 

V(/>f(so,si) = V[Li](so,si) 

Quadrilateral Eaee 

4>?j{so,si,to,ti) = (l)f{so,Sl)4>f{to,tl) 

= (l)f{so,si)V(j)f{to,ti) + (l)f{to,ti)V(l)f{so,si) 

Triangle Eaee 

4’ij{s0, Si, S 2 ) = </>f (so, Si)[-L|*](sO + Si, S 2 ) 

Si, S 2 ) = [L^*](so+si, S 2 )V(j)f{so, Si) + (/>f (so, si)V[L^*](so+si, S 2 ) 

A = l,2,3 

i = 2,...,ps 

iV = 2,3 

i = 2,...,ps 

j = 2,...,Pt 

iV = 2,3 

i>2,j>l 

i+j = 3,...,ps 

H (curl)* 

Operator 

Indiees 

Edge^ 

-Ef(so,si) = [Pi](so,si)(soVsi - SiVso) 

VxE;f(so,si) = (i + 2)[Pi](so,si)VsoxVsi 

Quadrilateral Eaee 

E;R(so,si,fo,ii) = </>f(to,ti)-Ef (so,si) 

VxE;R(so,si,fo,ii) = (ff{to,ti)VxEf{so,si) + V(/>^(to, ii) xElf (sq, si) 
Triangle Eaee 

E^{so,si,S2) = [L2^^](so+si,S2)Elf(so,si) 

VxE;,^(so,Si,S 2) = [L|^^](so + Si,S 2 )Vx£'f(so,Si) 

+V[Lj^^](so+si, S 2 ) X Ef{sQ, si) 

iV = 2,3 

1 = 0,... ,Ps-l 

N = 2,3 

1 = 0,... ,Ps-l 

j = ‘2,...,Pt 

N = 2,3 

i+j = l,...,Ps-l 

* For A'" = 2 the curl and cross product are V x ^ j and ( b) ) x ^ ^ = E 1 F 2 — 

0 In some cases the curl vanishes as shown in (4.10). 

E 2 F 1 respectively. 
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Ancillary Operators 

H{div) 

Operator 

Indices 

Quadrilateral Face^^ 

N = 3 

= Ef{so,si)xEf{to,ti) 

1 = 0, ... ,Ps-l 

V-V^{so,si,to,ti) = Ef{to,ti) • (VxE'f (so,si)) 

j = 0,...,pt-l 

- £’f(so,si) • iVxEf{to,ti)) 


Triangle Face^ 

N = 3 

Fif (so,si,S2) = [Pi]{so,si)[Pf"^^]{so+si,S2)[soVsixVs2 

i>0,j>0 

+ SlVs2 X Vso + ■S 2 VS 0 X Vsi^ 

i+j = 0,...,Ps-l 

V-Fif (so,Si,S 2 ) = ii+j + ^)[Pi]{so,si)[Pf''^^]iso + si,S 2 )Vso- (VsixVs 2 ) 


f In some cases the divergence vanishes as shown in (6.14). 

1: In some cases the divergence vanishes as shown in (7.14). 
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E.4 Quadrilateral 


Quadrilateral 


Geometry 


6 


4 


3 


1 2 
Affine Coordinates 

= 1-6 

V4^ = (-') 

/^o" = 1-6 

V/io^^ = (A) 



The order in the pair = {fx ^^, ) is pi. 

The order in the pair = (p,g^, pp) is p2- 




Shape Functions 

Indices 

Vertices 



a = 0,1, 6 = 0,1 

Edges 


<Pi = 

v<^i = (4“) + (4nv//|‘ 

L . . . ^ Pa 

(a, 6) = (1,2), (2,1) 

c = 0,l 

Face 


= 6y(/lo6/loi) 

* = 2,... ,pi 

v4 = v6°(4646 

j = 2,...,P2 
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Edges 
Et = 

VxEt = VntxEf{fl^l) 


Quadrilateral 


H (curl)"* 


Shape Functions 


Indices 

1 = 0, . . . ,_Pa-l 
(o, 5) = (1,2), (2,1) 
c = 0,l 


Family I 

4 = i?°(4\4^) 
vx4 = vx4(4,4) 

Family II 


4 = 

vx4 = vx4(4,4) 


H{div) 

Shape Functions 

For any given topological entity, the H (div) shape functions are the rotation of 
the corresponding FT (curl) shape functions: 


z = 0,... ,pi-l 
j = 2,...,P2 

i = 0,... ,P2-1 

j = 2,...,pi 


Indices 


V-F = VxE 


Shape Functions 


Indices 


i’lj = lfll(^o'i)WI(4’)(V«f‘ X V;.,= 


z = 0,... ,pi-l 

j=0,...,P2-l 


* In 2D the curl and cross product are V x ^ — §|^ and f ^2 ) ^ ( F 2 ) = E 1 F 2 — E 2 F 1 respectively. 




E.5 Triangle 


Triangle 


Geometry 



Affine Coordinates 

1^0 = I - Xi - X2 Ui = Xi 1/2 = X2 


The order in the niplet z7oi2 = {vq, i^i, 1 ^ 2 ) and all its subpairs is p. 




Shape Functions 

Indices 

Vertices 


< 

II 

a = 0,l,2 

= VVa 


Edges 



i = 2,...,p 

V<^| = V-^f (E,b) 

0 < a < b < 2 

Face 



* > 2, j > 1 

= V</.g(Eoi2) 

i+j = 3,...,p 
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Triangle 


Edges 

Ef = EfilVab) 
VxEf = VxEf{i^ab) 


H (curl)"* 


Shape Functions 


Indices 


0 < a < 6 < 2 


Family I 

4 = Et,{Voi^) 
VxEfj = VxEf^{i7oi2) 
Family II 

Eij = Efj{ivi2o) 
VxElj = VxEf^{i7i2o) 


H{div) 


Shape Functions 

For any given topological entity, the H (div) shape functions are the rotation of 
the corresponding FT (curl) shape functions: 


* > 0, j > 1 


Indices 


V-F = VxE 


Shape Functions 


i^ij = + Vl, Z/2)(Vl/i X Vj^2) 


Indices 

i > 0, j > 0 

z+j = 0, 


* In 2D the curl and cross product are V x ^ and f ^2 ) ^ ( F 2 ) = E 1 F 2 — E 2 F 1 respectively. 


■ 

I 


B] 




E.6 Hexahedron 


Hexahedron 


Geometry 



The order in the pair = (/Tq ^) is pi. 
The order in the pair /2 q^ = (p,g^, pp) is p 2 - 
The order in the pair /2 q^ = (p,g^, is pg. 
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Hexahedron 


-ff(curl) 


Shape Functions 


Edges 


VxEt= (4v4 +4V;u|^) xi^f(4) 

Faces 
Family I 

4 “ (^01 ’ 4) 

V X E^j = F-J'' V X Ej^j , /Igl') + X (/Igj", /Igl) 


Family II 


4 “ (4 ’ Foi) 


V X Elj = 4v X EYMm^fl^l) + x EYM^l^fl^l) 


-e 




Interior 
Family I 




^ ^ Fjfc 

Family II 




pd/j-fil ,-7'?2 


E/',-t?3N „ i7.n^-7-?i ,-r6^ 




V X £g» ^ 
Family III 




pO/-;^2 ,-7>?3 


NX fEI /-76 --€3^ 


= </>fc (4) 4 (4 > 4) 


Vx4, 




rat,7^3 ,-;?i 


E77?2\ I?n77?3 ,7?!' 


Indices 

1 = 0, . . . ,Pa-l 
(a, 6, c) = (l,2,3), 
(2,3,1), (3,1,2) 
d = 0, 1, e = 0,1, 


i = 0,... ,Pa-l 

j = 2,...,pb 

{a,b, c) = (l,2,3), 

(2, 3,1), (3,1,2) 
d = 0,1 

i = 0,... ,P6-1 

j = 2,...,pa 

{a,b, c) = (l,2,3), 

(2, 3,1), (3,1,2) 
d = 0,1 

z = 0,... ,pi-l 
j = 2,...,P2 
k = 2,...,p3 

i = 0,... ,P2-1 
j = 2,...,P3 
k = 2, ...,pi 

z = 0,... ,P3-1 
J = 2,...,pi 
k = 2, ...,P2 
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^ij ~ (/^Ol > /^Ol) 

Interior 
Family I 

— '/’f (/^Ol ) (/^oi ’ /^oi) 

Family II 

Family III 

y};k = ^l{4i)yIji4hfo\) 
v-F,5, = v</)f(4f)•F°(4^4^) 

L2 


Shape Funetions 


Interior 

= [«l(4‘)[-fil(Mo;)[ftl(4f)(Vpf‘ X 


Indiees 

1 = 0, . . . ,Pa-l 
j=0, ...,Pb-l 
(a,6,c) = (l,2,3) 
(2,3,1), (3,1,2) 
d = 0,1 

z = 0,... ,pi-l 

j = 0,...,p2-l 

k = 2,...,p3 

i = 0,... ,P2-1 
j = 0, ...,P3-1 

k = 2, 

i = 0,... ,P3-1 
j = 0,...,pi-l 
k = 2, ...,P2 


Indiees 

z = 0,... ,pi-l 

j=0,...,P2-l 
k = 0 ,... ,»3-l 




E.7 Tetrahedron 


Tetrahedron 


Geometry 


Ao = 1 - Xi - X 2 - X 3 
VAo= f-1 



Ai = xi 

VAi = f 0 

Vo 


A2 = X2 

VA2 = 

Vo 


As = X 3 


VAs = ( 0 


The order in the quadraple A0123 = (Aq, Ai, A2, A3) and all its subtuples is p. 
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Tetrahedron 

H (curl) 

Shape Functions 

Indices 

Edges 


Et = EfiXab) 

l = 0,...,p-l 

VxEt = Vx EfiXab) 

0 < a < 6 < 3 

Faces 


Family I 

1 > 0, j > 1 

4 = Ef^iXabc) 

z+j = l,...,p-l 

VxEIj = VxE^{Xabc) 

0<a<b<c<3 

Family II 

i>0,j>l 

Elj = Ef^iXbca) 

i+j = l,...,p-l 

VxElj = VxEf^{Xbca) 

0<a<b<c<3 

Interior 


Family I 


4» = [Lf+">](l-A3.A3)£$(X„i2) 

z > 0, j > 1, fc > 1 

= [Lf+^)](1-A3, A3)VxE,^(Aoi2)+V[Lf+^'^](l-A3,A3)xF;,^(Aoi2) 

i+j+k = 2,...,p-l 

Family II 



1 > 0, j > 1, fc > 1 

VxF;^, = [Lf+^)](l-Ao, Ao)VxE,^(Ai23)+V[Lf+^'^](l-Ao, Aq) xF;,^(Ai23) 

i+j+k = 2,...,p-l 

Family III 


^*yfc = [^f^'^](l-Ai,Ai)i?§(A23o) 

i>0,j>l,k>l 

VxF;,^., = [Lf+^)](1-Ai, Ai)VxE,5(A23o)+V[Lf+^'^](1-Ai, Ai) xF;,^(A23o) 

i+j+k = 2,...,p-l 
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Tetrahedron 


H{div) 


Shape Functions 

Indices 

Faces 

i > 0, j > 0 

Vlj = V^iXabc) 

i+j = 0, 

y-v!^ = y-v^iKUx)) 

0<a<b<c<3 

Interior 


Family I 



i>0, j>0, k>l 

V-rAt = [Lj“^‘>|{l-A3,A3)V.V'f(A„,2) + V|4l't^'>](l-A3,A3)-V'if(X„i2) 

i+j+k = l,... ,p-l 

Family II 


Kf = [ii'"^'’l(l-Ao, A„)FA(Ai23) 

z > 0, j > 0, fc > 1 

V-V'5t = [i.?‘^‘’l{l-Ao,A(,)V.V'A(Ai2s)+V|4''k^'>](l-A„,A„).VAf(X,23) 

z+j+fc = l,... ,p-l 

Family III 


Kf = [ii'"^'’l(l-Ai, Ai)FA(A23„) 

z > 0, j > 0, fc > 1 

V-r.5t = [L^“^‘']{l-Ai,Ai)V.\Af(A2s„) + V|4l't«-'>](l-A„Ai).V'f(X23o) 

z+j+fc = l,... ,p-l 


L2 

Shape Functions 

Indices 

Interior 

= [fll(Aoi)|Z’2"*‘l(Ao+A, 

, X2)[p!^^\i-X3, A3)(VAi X VA2)-VA3 

i > 0, j > 0, fc > 0 

i+j+k = 0,... ,p-l 
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E.8 Prism 


Prism 


Geometry 



Uq = 1 - Xi - X2 

Vl = Xl 

l'2 

/ -1 \ 

f l\ 


< 

II 

1 

Vvi = ( 0 

Vl'2 

Vo/ 

\oJ 



^io = l- z 
V/Uo = ( 


^Jix= z 
V|Ui = 


= 1 


The order in the triplet r'oi 2 = (r'o, ^^ 2 ) and all its subpairs is p. 

The order in the pair /loi = (Mo: Mi) is 9- 
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Prism 

H (curl) 

Shape Functions 

Indices 

Interior 


Family I 

Eijk = <Pkim)Efj{uoi2) 

* > 0, j > 1 

= 4>^{floi)VxEf^{i7oi2) + V((>f (^oi) x£;,^(Foi2) 

k = 2, ...,q 

Family II 

i>0,j>l 
i+j = l,...,p-l 

Eijk = </>fc (/^oi)-E'jy (z^i2o) 

= (p^{iloi)V X Efj{Pi 2 o) + V(()f (/lot) x£')^(Fi 2 o) 

k = 2, ...,q 

Family III 

i > 2, j > 1 
i+j = 3,...,p 

Eijk = (moi) 

V X E^jf^ = V(/>g(Foi 2 ) X E^ipoi) 

k = 0,..., q—1 
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Prism 

H{div) 

Shape Functions 

Indices 

Triangle Faces 

(F 012 ) 


i > 0, j > 0 
i+j = 0, 

V-F,f =V/re-F,f(Foi2) 


c=0,l 

Quadrilateral Faces 

^ij = V^{^ab,Poi) 


II II 

0 0 

1 1 

V-V!j=V-V^{Uab,m) 


0 < a < 6 < 2 

Interior 



Family I 


* > 0, j > 1 
z+j = l,...,p- 



V-F,^, = Ef(Moi)-(vxi?§(Foi2)) 


k = 0,..., q—1 

Family II 


i>0,j>l 

^ijk = (z^i2o) X-E'fc (Foi) 


V-F,^, = Ef(Moi)-(vxi^,5(Fi2o)) 


1 

cT 

II 

Family III 


i > 0 , j > 0 

z+j = 0, 

^^jk = ( 1 ^ 012 ) 


V-F,^, = V,/.f(/2oi)-F,f(Foi2) 


k = 2, ...,q 

L2 

Shape Functions 

Indices 

Interior 


i > 0, j > 0 

i^tjk = [Pi]im)[Pf'^^]{i^o + i^i,i^2)[Pk]im)iVi^ixViv2)-Vfii 

z+j = 0, 




E.9 Pyramid 


Pyramid 


Geometry 



= ( -1 



Affine Coordinates 


= 1 - 


fi 

i-C 

V/^o’ ^ = n-n2 ( 0 

V -6 

$2 
1-C 


, C:?l _ 

0 “ (T^ 

= 1 - 


\7,, _ 1 

V / Tq — ( 1_^)2 

/^O = 1 - C 


V/Tq — 





,C-0 - 


= 6 




, C’O — 


= 6 


Vz/f = 


0 

-(1-C) 

-C2 


,C-0 — 


= c 


= 


,C-0 _ 


= c 


Vz2, 


C,C2 _ 


= ,s 


(1-C) 

0 

Cl 


,,C:Cl _ Cl 

Ml — l_f 

V7,,C:Cl _ 1 

VMi - (ir^ 

,,C:C2 _ C 2 
Ml — 1-f 

VmP = ,t^((A)) 

Ml = C 

V/i« = 


Pyramid “Affine” Coordinates 


Ai = 


_ (1-Ci-C)(i-C2-C) 


1-C 


A 2 — 


_ Ci(i-C2-C) 


1-C 


A 4 = (^-^/_-0C2 


A5 — c 


A3 = !$ 



/-(i-C2-C)\ 


/ (1-C2-C)\ 


( C2 \ 

VAi = 

1-C 

-(I-Ci-C) 

1-C 

C 1 C 2 1 

\(i-c)2 V 

VA 2 = 

1-C 

-Cl 

1-C 

-C 1 C 2 

v (i-c )2 y 

VA3 = 

1-C 

Cl 

1-C 

C 1 C 2 

\(1-C)V 


VA 4 = 


/ ^ \ 

1-C 

(1-Ci-C) 

1-C 

-C 1 C 2 

V (i-c)2 y 


VAs = 


Col 

0 

vv 


The order in aff tupfes of coordinates is p. 
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Pyramid 


Shape Functions 

Indices 

Vertices 


r = K 

0 = 1 , 2 ,3,4,5 

V</>" = VAa 

Mixed Edges 

i = 2,...,p 

( 0 , 5 ) = (1,2), (2,1) 

V4>t = * 

c=0,l 

Triangle Edges 

</>! = (pfiKb) 

V<PI = VcPfiXa^) 

Quadrilateral Eace 

i = 2,...,p 

a = 1,2, 3,4 

4>lj = Pmh 

i = 2,...,p 

V4 = + </>°(/2o¥\/2o¥^)Vmo^ 

j = 2,...,p 

Triangle Eaces 

* > 2, j > 1 

4>lj = 

i+j = 3, ...,p 

( 0 , 5 ) = ( 1 , 2 ), ( 2 , 1 ) 

c = 0 ,l 

Interior 


4fjk = 4'?j{Poi' ^ Pkh4>k(.Pm) 

1 — 2^... 

j = 2,...,p 

V</>4 = 4'k (Foi ) ^^?j (Fo4' > ) + ^?j (Foi ) 

c^r 

II 
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Pyramid 

H (curl) 

Shape Functions 

Mixed Edges 

Et = 

VxEt = V X Efii7,f^) + X (4«“) 

Triangle Edges 

Et = EfiXa,) 

VxEt = VxEf{Xa5) 

Quadrilateral Eace 
Family I 


4 = 


Vx4 = + 

Family II 


,-rC.6 


Cv7,,C s. pn/'.-tC:?! ,-fCI?2^ 


4 = (Ml)"4(4i'^”.4^‘) 


Vx4 = (4)VxE“(^«>,^«') + 2;,„^V,.ixE“(^«>,^«') 
Triangle Faces 
Family I 

vx4 = 4’«*>vxi?,44t) + x4(4t) 

Family II 

4 = 

vx4 = 4’«*>vx4(4«“) + v//^ x4(4«“) 




Cy7,,C^ IT’D/',-?CiS2 


Indices 
i = 0,...,p-l 

(a, 5) = (1,2), (2,1) 

c=0,l 

i = 0,...,p-l 
a = 1,2, 3,4 


i = 0,...,p-l 
i = 2, ...,p 

i = 0,...,p-l 

j = 2,...,p 


* > 0, j > 1 

z+j = l,...,p-l 
(a, 5) = (1,2), (2,1) 
c=0,l 

* > 0, j > 1 

z+j = l,...,p-l 
(a, 5) = (1,2), (2,1) 
c=0,l 
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Pyramid 

H (curl) 

Shape Functions 

Indices 

Interior 





Family I 



7—9 

..,p 

..,p 



t — z,. 

J = 2,. 

VxEfji^ = 0 



fc = 2,. 

■ ■ ,P 

Family II 



1 = 0,... 

J = 2,. 

.,p-l 

..,p 




VxE^, = (4)VxS°(4i«S4«^) 



k = 2,. 

■ ■ ,P 

+ 4>k iPoi) + 4>k iPm)^f^o) 

1 x4(4^'^ 




Family III 



1 = 0,... 

J = 2,. 

,,p-l 

..,p 

E^k = 



vx4, = (4)vx4(4«^4«^) 



k = 2,. 

■ ■ ,P 

+ (^l^o'^^k (4) E 4>k (Foi)^Fo^ 

1 x4(4^"^ 

rof ) 



Family IV 



A _O 

..,p 

..,p 




Z — ^ 7 * 

J = 2,. 

VxEj; = 



n = max{l, j} 
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Pyramid 

H{div) 



* To avoid computing (//f^ and see (B.42)-(B.45) for alternate expressions. 
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Pyramid 

H{div) 






E.IO Orientations 


Orientations 


Local Orientations 



a 



Triangle 




Hexahedron Tetrahedron 




Pyramid 
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Edge 

(«0,Sl) = 

Quadrilateral Face 


Triangle Face 


O'o (so, Si, S2) = 


Orientations 

Focal-to-Global Permutation Functions 


Clo (so,Sl,to,G) = 


o-f(so,si) = (so. Si) 
O'f(so,si) = (si,So) 

(T^(so,Si,to,G) = (s 
(t 5 ^(so, Si,to,G) = {t 
0-2 (so. Si, to, G) = (s 
0-3 (■So, Si, to, tl) = (t 
0-4 (so, Si,to,ti) = (t 
0-5 (so. Si, to, ti) = (s 
0-6 (so. Si, to, ti) = (t 
(T^(so, Si,to,ti) = (s 


if o = 0 
if o = 1 


(so, 

Si, to 

^l) 


if 0 = 

= 0 

(fo. 

tl, Si, 

So) 


if 0 = 

= 1 

(si, 

So, tl 

to) 


if 0 = 

= 2 

(fl. 

to. So, 

Si) 


if 0 = 

= 3 

(fo. 

tl. So, 

Si) 


if 0 = 

= 4 

(si, 

0 

0 

^1) 


if 0 = 

= 5 

(fl. 

to, Si, 

So) 


if 0 = 

= 6 

(so, 

Si, tl 

to) 


if 0 = 

= 7 

), Si 

S 2 ) 

if 

0 

= 0 


L, S2 

so) 

if 

0 

= 1 


J, So 

Si) 

if 

0 

= 2 


), S2 

Si) 

if 

0 

= 3 


L, So 

S 2 ) 

if 

0 

= 4 


>, Si 

So) 

if 

0 

= 5 



In the expressions for the shape functions, precompose the ancillary operators (and their 
differential form) with the corresponding permutation function to obtain orientation em¬ 
bedded shape functions. 






